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FOREWORD 


The  TRIOIL  computer  code  described  herein  is  as  it  existed 
on  July  1,  1967.  The  code  has  been  in  continuous  development  for 
three  years  and  in  its  presented  form  has  been  applied  in  this  re¬ 
port.  However,  the  development  and  improvement  in  both  the  physics 
and  mathematics  of  the  code  are  being  continued,  so  that  duplication 
of  results  (or  even  close  agreement)  between  problems  run  with  the 
code  as  published  and  the  code  as  it  existed  either  before  or  after 
this  time  is  not  necessarily  to  be  expected. 

General  Atomic  has  exercised  due  care  in  preparation,  but 
does  not  warrant  the  merchantability,  accuracy,  and  completeness 
of  the  code  or  of  its  description  contained  herein.  The  complexity 
of  this  kind  of  program  precludes  any  guarantee  to  that  effect. 
Therefore,  any  user  must  make  his  own  determination  of  the  suita¬ 
bility  of  the  code  for  any  specific  use  and  of  the  validity  of  the 
information  produced  by  use  of  the  cple. 


U^unpnt  ha,,  ^ 
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1.  INTRODUCTION 

The  TRIOIL  code  is  a  three-dimensional  continuous  Eulerj.an  cartesian 

*  i 

hydrodynamic  code.'  It  is  a  natural  extension  of  the  two-dimensional  OIL 

code  (Ref.  l).  The  code  is  in  an  infant  stage  as  of  now  \>ut  work  is  con¬ 
tinuing  to  make  it  an  operative  tool  for  research  in  the  field  of  msnerical 
hydrodynamics.  In  this  version  of  TRIOIL,  the  scalar  pressure  is  the  only 
stress.  However,  it  appears  to  be  straightforward  to  incorporate  stresses 
due  to  strength  and  viscous  forces  as  has  been  done  for  the  two-dimensional 
version  (Ref.  2)  and  adding  two  materials  (Ref.  4). 

The  maximum  size  of  the  three-dimensional  grid  is  limited  by  the 
present  32K  to  64K  storage  computers.  In  this  present  version,  the  maximum 
number  of  cells  in  any  single  direction  *IMAX  for  the  x-direction,  JMAX  for 
the  y-direction  and  KMAX  for  the  z-direction)  is  limited  to  30>  further,  the 
total  number  of  cells  is  limited  to  6000.  But,  with  high  speed  disks  and 
drums,  as  are  currently  being  made  available,  one  can  store  x-y  slabs  of 
data.  This  means,  at  any  given  time,  we  have  in  core  memory  only  two  slabs 
of  x-y  data,  which  is  sufficient  for  the  various  phases  of  the  program. 

This  procedure  will  essentially  remove  the  limitations  on  grid  size,  and 
enable  one  to  apply  this  code  to  solving  many  types  of  three-dimensional 
problems. 

A  rezone  routine,  that  increases  all  linear  cell  dimensions  by  a 
factor  of  two,  is  presently  an  operative  feature  of  this  code.  This  routine 
is  programmed  primarily  for  hypervelocity  impact  calculations  into  an  in¬ 
finite  target.  To  use  it  for  different  applications  will  require  some 
minor  modifications.  Eight  ceils  in  the  old  grid  are  combined  into  one  for 
the  new  grid.  Energy,  mass,  and  momentum  are  conserved.  This  routine  as¬ 
sumes  the  following: 

(a)  That  IMAX,  JMAX  and  KMAX  are  even  integers.  (iMAX,  JMAX 
and  kMAX  are  the  total  number  of  cells  in  the  x,  y,  and 
z  directions. ) 

* 

This  work  was  funded  Jointly  by  the  Systems  and  Research  Department 
of  Honeywell,  Inc.,  at  Minneapolis,  Minnesota  and  General  Atomic  Division 
of  General  Dynamics  Corporation. 
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(b)  That  the  J  value  of  the  initial  interface  between  pro¬ 
jectile  and  target  is  2^,  allowing  N  rezones. 

This  routine  can  be  called  for  whenever  mass  leaves  any  of  the  six 
grid  boundaries. 

The  six  slides  of  the  grid  can  either  be  reflective  or  transmlttive. 
Variable  zoning  of  Ax,  by  and  hz  is  a  working  feature  of  the  code. 

Active  grid  counters  are  calculated  for  all  three  directions.  By  doing 
this,  one  processes  only  that  portion  of  the  grid  that  is  active. 

Programs  for  displaying  density,  velocities  and  pressures  are  cur¬ 
rently  being  programmed  using  the  Stromberg-Carlson  4020  Plotter. 

A  subroutine  called  SETUP  is  available  for  generating  the  starting 
data  for  the  TRIOIL  code.  This  routine  will  only  generate  a  rectangular 
parallelpiped  for  the  projectile  and  target.  For  the  more  complicated 
geometries,  we  plan  to  modify  the  CLAM  code  (Ref.  l)  to  generate  three- 
dimensional  starting  conditions  for  the  TRIOIL  code. 

An  accurate  estimate  of  computer  time  required  to  run  a  realistic 
problem  in  three  dimensions  is  not  presently  available,  however,  calculations 
reported  here,  have  been  run  on  the  UNIVAC-1108  and  Control  Data  6600 
computers  to  a  point  where  the  peak  pressures  have  been  attenuated  by  se¬ 
veral  orders  of  magnitude  for  computer  times  of  approximately  one  hour. 


T2_e  code  is  written  in  FORTRAN  IV  language  and  ic  operative  on  the 
IBM-7044  and  the  UNIVAC-1108  as  well  as  the  Control  Data  6600. 

2.1  Basic  Equations 

The  Eulerian  equations  in  cartesian  geometry  are: 

v  '  at  ax  ay  az 

(B)  P  +  pu  §»  +  pv  |»  +  Pw  |J  +  g  -  0  . 


^ $ +  »'u !; * (,v  w *  o'* +  w*° 


V 
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dw  dw  3w 

p  57  +  pu  57  +  pv  5y  T  •*-  &z  ’  az  -  ' 

SB  5E  3E  5E  5(Pu) 

P  St  +  pu  57  +  pv  57  +  pw  5T  +  “Tx 


<c)  p  If 


These  equations  are  solved  In  two  parts,  as  in  the  OIL  code  (Ref.  1) 
and  the  familiar  particle-in-cell  codes  (Ref. .3)*  The  transport  terms: 
on  the  left  side  of  Eq.  (B)  and  (C)  are  temporarily  dropped,  while  we  com¬ 
pute  (in  PHI)  the  ipomentum  and  energy  contributions  due  to  pressure  forces 
only.  The  omitte^  transport  contributions  to  the  momentum  and  energy  are 
later  approximated  when  we  solve  equation  (A)  and  move  mass  across  the  cell 
boundaries. 


2.1.1.  Effect  of  Pressure  Forces  Only  (PHI)  . 

Re-writing  Eqs.  (B)  and  (C)  with  the  transport  terns  dropped 


5u 


ap 


p  St  ■  •  55 


Sv  SP  \ 

p  5t  "  “  5y 

_  aw  ap 

p  at  ’  ’  57 

aE  apu  apv  apw 

p5t  "  _  57”  "  57"  ’  5T" 

p  *  f(p,l)  equation  of  state 
These  variables,  in  one  consistent  set  of  units  are 


p  ■  density  of  cell  (L)  in  g/cm^. 
x  ■  x  coordinate  in  cm. 
y  ■  y  coordinate  in  cm. 
z  -  z  coordinate  in  cm. 
u  ■  x  component  of  velocity  in  cm/ shake, 
v  -  y  component  of  velocity  in  cm/ shake, 
w  ■  z  component  of  velocity  in  cm/ shake. 


(1) 

(2) 

(3) 

(»0 

(5) 
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P  ■  material  pressure  in  Jerks/cm^, 

E  m  total  specific  energy  in  Jerks/g. 

I  ■  specific  internal  energy  in  Jerks/g 
t  ■  time  in  shakes 
1  shake  ■  10~®  sec 
1  Jerk  ■  10A^  ergs 

There  are  no  built-in  units  for  this  code,  the  user  can  select  his 
own  units  in  a  given  application  by  way  of  the  equation  of  state  constants 
and  input  data. 

The  density,  pressure,  velocities  and  internal  energy  are  all  cell- 
centered  quantities  referred  to  with  the  index  (L). 


Rewriting  Eq.  4: 


3E  dPu  dPv  3Pw 

p  St "  "  5x"  *  5F  '  5T 


a  rT  i/  2  2  2vi  apu  apv  apw 
p  5t  l1  +  *<u  +v  +w  >J  "  -  55T  -  ST  ’  5T 


ai  au  Sv  ^  aw  _  au  „  ap 

p57  +  pu5t  +  pv5t  pwSt“"p5x"  5x 


au  ap  p  av  ap 
5x*“5x'PSy  3y 


p|s.  „j£ 

*  az  OZ 


thus 


„  &U  &P  „  dV  3P  «  SW  &P 
p  5t  “  "  Sx  and  p  St  "  -  5y  and  p  St  *  "  H 


ai  _  r au  av  aw 
p  8t  ■  •  p  [e  *  57  +  Si 


Now  we  can  write  the  four  differential  equations  (1,  2,  3  and  6)  in 
difference  form. 

The  storage  arrays  for  the  cell  centered  quantities  (mass,  velocities, 
pressure  and  specific  internal  energy  are  as  follows  (Fig.  1). 
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In  the  discussion  to  follcw,  please  refer  to  Fig.  2.  Re-writing  the 
equations  (l  to  3)  in  difference  form  results  in: 

The  x-momentum  equation  (l)  becomes  in  difference  form 

IvL  -  U?\  Pn  -  Pn 

n  j_L _ L-l  L+l 

PL  l  4t  I"  2Ax(i) 

The  y-momentum  equation  (2)  becomes  in  difference  form 

n 
pL 

The  z- momentum  equation  (3)  becomes  in  difference  form 

(w  n\  pn  pn 

L  "  L 1  L-(iMAX)(.1MAX)  '  rL+(lMAX)( JMAX) 

At  j  "  2Az(k) 

Here  the  acceleration  of  cell  L  is  seen  to  depend  only  on  the  pressures  in 

the  neighbor  cells  (not  that  of  L)t  Defining  pressures  at  interfaces, 


At 


pn  -  pn 

L-iMAX  L+iMAX 

2Ay(J) 


PLn 

PRRn 

PBLOn 

PABOVEn 

PBINDn 


p*1  4.  pn 

PL-1  *L 

2. 


Pn  +  Pn 
L  rL+l 


2. 


P11  +  p° 

L-iMAX  rL 

2. 


Pn  +  Pn 
L  L+1MAX 

2. 


L-(iMAX)(.1MAX)  L 

2. 
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and 


VBLO  - 


v  +  v 
L  L-iMAX 

2. 


VABOVE 


v  +  v 
L  L+iMAX 

2. 


UL 


UL  +  UL-1 


URR  = 


UL  +  UL+1 

2. 


UBIND  «=  WL  *  WL-(1MAX)(JMAX) 
2  • 


WZR 


WL  *  WL+(iMAX)(.1MAX) 

2. 


then 


n 

pL 


17 '  -  i?  I 

L  L 

=  Pn 

At 

M 

VBLOn  VbEo 
2Ay(j) 


VABOVE*1  +  VABOVE  ULn  +  UL 
2Ay(j)  2AX(i) 


URRn  +  UKR  .  UBIND n  +  UBTND  WZRn  +  WZr] 

2AX(i)  2AZ(k)  2AZ(k)  J 

The  solution  of  the  momentum  equations  arc  very  straightforward,  how¬ 
ever,  the  solution  to  the  energy  equation  requires  the  velocities  at  two 
different  time  steps  (the  old  and  new  velocities).  We  have  chosen  to  make 
two  passes  through  the  routine,  the  first  pass  to  integrate  the  momentum 
equations,  but  formulate  the  interface  velocities  first  for  the  work  term 
contribution  to  the  internal  energy  before  integrating  the  momentum  equations 
(since  we  have  allowed  only  one  array  per  velocity  component). 

The  second  pass,  we  bypass  the  momentum  equations,  and  only  compute 
the  new  interface  velocities  for  their  contribution  to  the  work  term.  A 
single  pass  can  be  done  by  looking  ahead  two  cells  above;  two  cells  to  the 
right  and  two  cells  in  front. 
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The  choice  of  the  velocities  at  n  +  £  in  the  energy  equation  is 
apparent  in  the  following  discussion.  For  convenience  ve  will  go  through 
the  logic  in  the  y  direction.  Since  we  have  dropped  the  transport  terns, 
our  integration  of  the  momentum  and  energy  equations  have  not  been  advanced 
to  time  (n  +  l).  As  before,  we  designate  the  phase  1  velocities  and  energy 
as  ~,  w,  and  If. 


and 


+ 


+ 


At 

W-3/a  - 

n 

PH 

- 

2Ay(J)  J 

itprt 

vJ-3/a  '  L*l 

pa 

L 

j 

where 


VJ-3/2 


V3/2  *  VJ-3/2 

2. 


and 


2. 


Before  entering  PHI,  where  the  quantities  are  at  time  by  the  totfil 
energy  of  the  system  (considering  one  dimension,  the  y  direction)  is: 

JMAX 


'  *  t  [lU +  *  (vh]1  • 


And  the  total  energy  at  the  end  of  PHI  is  then 


The  change, 

AE  » 


AE  =  E11  -  If  should  be  zero  for  energy  conservation  or 
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the  A  kinetic  terms  can  be  represented  by 

lvIU ;  ial 


2. 


JMAX 

“■  z 

J=1 


"“j-i  [tj-4  -  Vi +  vj-i  (v"-i  '  Vi> 


JMAX 

I 

J«1 


AMX 


j-i 


At  1  j-k 

V3/2  -  Vfel 

n 

l  2  Ay  1 

L  pJ-i 

J 

-  V 


£1_  n-3/a  :  PJ<i 


*>U\  2 ^ 

***  AMXn 


At 


ZflflA.  i  r 

- -  ^  —  I  ■ -  pn  1  V  / 

j-l  0 u  2*yJ  1  *  J'3/2 


+  PJ-i  vJ*i  '  *^-3  vJ-i  +  PJ-*i  vJ-iJ 


JMAX 


DX(  i)DZ(  k)  £ 
J-l 


J4  VJ-3/2  +  j-3/2  VJ-i 


VJ-i  *  PJ-£  VJ^ 


Thu6,  the  last  two  terms  in  J  being  cancelled  by  the  first  two 
terms  in  J+l.  By  prescribing  the  proper  boundary  conditions,  we  will  have 
energy  conservation  for  the  entire  grid. 


EXAMPLE: 


for  J  ■  1 


Note,  the  boundary  terms  do  not  cancel.  For  the  example,  assume  the 
bottom  boundary  to  be  reflective.  Referring  to  Eq.  (A),  the  first  two  terms 

■4 


do  not  cancel.  If,  however  we  set  the  pressure  of  the  mirror  cell  P 


,n 

■i 


(this  does  not  imply  that  v^  =  0)  and  v_^  =  -y^, 
lation  of  the  first  two  terms.  If,  however,  we  designate  the  bottom  boundary 
to  be  transmittive,  our  boundary  conditions  are 
p"/„,  and  that 


thiB  does  lead  to  cancel- 


0.  which  implies 


Q  ‘  r3/2’ 


i  "  ”3 

This  still  leaves  us  with  the  two  terms 


Y  ~  Y 

-r  ^/2  *  thus  we  c0mPen8a'te  for  these  terms,  by  adding  or  sub¬ 

tracting  energy  to  the  syi  tern.  To  keep  all  the  books  straight,  E  th  must 
be  modified  accordingly.  Similar  boundary  conditions  may  exist  for  the 
other  four  sides  of  the  grid. 


2.1.2  Adding  the  Transport  Terns  (PH2) 

Rewriting  Eq.  (A)  the  mass  transport  equation  in  finite  difference 
form  results  in 


n+1 


“  Pi 


n 


At 


n  ~ 


n  ~ 


n  ~ 


n  ~ 


n  ~ 


PJ-1  Vi  •  P.1  V.1  t  pi-i  Ui-1  •  pi  ui  .  pk-i  Vi  -  Pk  w] 

■*—  ■  n  "■  ■  11  +  "  ■  ■  ■  1  "  + - *  -  ■— 


Ay(  J) 


AX 


U) 


AzTkT 
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or  In  the  formulation  for  mass  we  have 

AMl  =  At  £(Av)^p  -  (Av)^p  +  (Au)£p  (Au)^p  +  (AwJ^p  -  (Av)^pJ 

where  the  superscripts  X,  Y,  Z,  refer  to  the  particular  coordinate,  and 
B  =  bottom,  T  =  top,  L  =  left,  R  =  right,  BH  =  behind  and  F  °  front.  The 
mass  that  is  moved  from  cell  to  cell,  has  then  momentum  and  energy  associ¬ 
ated  with  it,  thus  these  are  the  approximations  to  the  transport  terms  that 
were  omitted  in  PHI  for  the  momentum  and  energy  equations. 

The  p  (density)  is  that  of  the  donor  cell,  and  the  velocity  is  yet 
to  be  determined.  Various  techniques  for  the  velocity  weighting  have  been 
tried  (Ref.  l).  The  velocity  weighting  scheme  in  this  report  is  identical 
to  that  in  the  OIL  report  (Ref.  l). 


2.1.3  Time  Control 

The  time  step  (At)  is  controlled  by  the  Courant  condition  and  the  con¬ 
dition  that  the  mass  flux  equation  will  not  empty  more  mass  from  a  cell  than 
is  there. 


Take  the  y  direction  (a  similar  treatment  is  done  for  the  other 
two  directions) 


Let 


then 


AM  =  pvAAt 

y 

"  "  V(L) 
p  *  P(L) 

AMy  -  AMX(l) 
AMX(L)  =  P(L)  V(L)  Aj 


HWA/.  x 

"  DX(i)  Dy|j)  DZ(k)  V(L)  DX^  132 ^  Lt 
^(L)  V(L) 

Dy ft) 
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0r  1v(l)I  ^  «))  such  that  the  flux  in  the  y  direction  will  not 

empty  the  cell  ( L) . 


The  Courant  condition  is  that  £  1  where  C  =  sound  speed  and  (L. 
is  the  maximum  speed  at  which  a  disturbance  can  propagate  in  the  given  grid 


or  Ct-<  1.,  or  1.,  and  C  7-  £  1. 

Ax  *  Ay  ’  Az 


2.1.4  Remarks 

No  corner  coupling  (that  is,  mass  is  constrained  to  move  at  right 
singles  to  the  sides  of  the  cell)  exists  in  this  version,  and  no  attempt  to 
systematically  study  this  has  been  initiated.  The  movement  of  mass  across 
the  cell  boundaries  give  rise  to  forces  which  are  effective  in  reducing 
fluctuations  that  arise  from  the  differencing  technique.  That  is  of  the 
form  of  a  "true"  viscosity,  being  proportional  to  the  velocity  gradient. 

It  is  tnis  force  which  enables  the  Eulerian  codes  to  treat  shocks,  where 
again,  as  in  the  Q  method  used  in  Lagrangian  codes,  the  shock  is  spread 
over  two  or  three  zones. 

2.2  Logic  of  TRIOIL 

The  logip  involved  in  following  a  given  cell  (L)  from  t  to  t  +  At 
or  cycle  n  to  n  +  1  is  as  follows: 

We  assume  we  have  integrated  the  mass,  the  three  velocity  components 
and  the  internal  energy  to  cycle  n,  now  all  that  remains  to  complete  cycle 
n,  and  to  begin  cycle  n  +  1,  is  to  update  the  pressures  from  the  equation 
of  state  and  calculate  a  new  time  step. 


2.2.1  CDT  Routine 

Here  we  calculate  the  pressure  (P)  array  for  the  entire  grid.  The 
pressure  ( P^)  =  f  (p^,  1^)  where  L  is  the  index  of  the  cell  in  question, 
defined  as  L  =  (j-l)  iMAX  +  i  +  (k-1)  (  iMAX)  ( JMAX) .  The  density 

several  times  during  the  cycle.  It  is  planned  in  the  future,  to  replace  the 

mass  storage  with  density.  The  speed  of  sound  C  *  for  a  polytropic 

ap  l/2  p 

equation  of  state)  or  (5—)  '  for  a  general  form  is  then  calculated.  From 

OP  8 

the  particle  velocities  and  speed  of  sound,  a  new  At  is  calculated. 


DX(i)  Dy^j)'  DZ(k)^  is  not  one  the  variab^-es>  i-fc  must  calculated 
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(Options  for  At  control  are  identical  to  that  in  OIL.)  The  cycle  number 
and  the  time  are  now  advanced. 

2.2.2  The  EDIT  Routine 

This  code  has  three  separate  editing  routines  all  included  in  the 
subroutine  EDIT.  The  first  of  these  is  the  routine  called  "Short  Print". 

This  displays  the  problem  number,  cycle  number,  time,  internal  and  kinetic 
energy,  energy  check  and  indices  of  the  cell  that  is  controlling  the  time 
step.  The  next  routine  available  is  called  "Long  Print".  Here  one  edits  each 
column,  the  three  velocity  components,  pressure,  mass,  density,  specific  in¬ 
ternal  energy  and  y,  all  versus  J  (the  index  of  the  rows).  Thus  there 
are  kMAX  sets  of  these,  beginning  with  k=l.  The  normal  units  chosen  are 
gram-cm-sh.  which  gives  a  logical  unit  for  the  pressure  as  10  ergs/cnr 
and  for  the  specific  internal  energy  as  10^  ergs/gram. 


2.2.3  PHI  Routine 


Here  we  integrate  the  three  momentum  equations  and  the  internal 
energy  equation  due  to  pressure  forces  only.  No  material  is  moved  at  this 
time,  and  the  transport  terms  are  dropped  temporarily.  Using  the  new  pres¬ 
sure  and  At,  we  can  solve  the  momentum  and  energy  equations. 


PL(j),  the  pressure  at  interface  (i-l)  and  UL(j),  the  velocity  at 
interface  ( i-l)  are  available  from  the  previous  column  sweep  on  i-l: 


PL(J) 


Pn  +  Pn 
L  L-l 


UL(J)  = 


^  +  u£-i 


The  FBLO  term  which  was  the  PABOVE  for  the  cell  below  (L-iMAX)  and 
VBLO  which  was  VABOVE  for  cell  below,  are  also  available  for  interface  J-l. 


PBLO 


pH  +  pn 

L  L-iMAX 


v“  +  V11 
L  L-iMAX 


VBLO 


2. 
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The  PBIND  term  (the  pressure  at  the  back  surface  of  the  cell)  which 
was  PZR  for  cell  L  -  (iMAX)  ( JMAX)  and  UBIND  (the  velocity  at  the  back  sur¬ 
face  of  the  cell)  which  was  WZR  for  cell  L  -  ( iMAX)  ( JMAX)  are  also  available 
for  interface  K-l. 


PBIND 


Pn  +  P° 

L  L  -  (iMAX)  (JMAX) 

2. 


UBIND  = 


L  -  (iMAX)  ( jMAX) 


Thus,  we  need  only  to  calculate  quantities  at  the  top,  the  right,  and 
the  front  of  cell  (L).  At  the  top  we  calculate 


PABOVE 


Pn  +  Pn 
L  L  +  iMAX 

2. 


v°  +  v°  pn  +  pn 

and  V ABOVE  «  — - ^  ;  at  the  right  we  calculate  PR  ■  — L 


and  URR  as 


2. 


2. 


PZR 


and  in  front  as 
Pn  .  pn 

rL  rL  +  (  iMAX)  (  JMAX) 

2. 


WZR  =  L  L  1  L  mAX)  (  JMAX) 
2  • 


When  the  cell  in  question  is  void,  the  pressures  at  the  top,  right 
and  front  interface  are  set  to  zero  and  the  velocities  are  set  =  to  the 
velocity  of  the  cell  above,  to  the  right,  and  in  front  respectively.  If 
a  occupied  cell  has  a  void  neighbor,  the  pressure  at  that  interface  is  set 
*  0,  and  the  velocity  at  that  interface  is  set  =  to  the  velocity  of  the 
occupied  cell  in  question. 

We  now  have  sufficient  information  to  integrate  the  three  momentum 
equations  and  part  of  the  internal  energy  equation. 


P 


Bu  -Bp 
Bt  "  Bx 


PL?  -  PKRn 

+  — -  Dy(d)  DZ(k)  At 


or 
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and 


or 


and 


P 


dv  _  -dp 

dt  "  dy 


V  =  V  + 

L  L 


PBLOn  -  P ABOVE11 


AMX 


n 


DX(  i)  DZ(k)  At 


P  Bt  dz 


~  -  (pBIND(M)n  -  PZgj^j)  ^  At 

AMx£ 


where  the  index  m  refers  to  the  ( i,  j)  value  of  the  slab  x*y  just  behind. 
We  can  add  the  work  terms  due  to  velocities  at  cycle  n  to  the  change  in 
internal  energy  as 


n  PL  At 
jn  +  k  _  .  _ 

L  Am,n 
AMX. 


(VBLOn  -  V ABOVE11) 

2. 


DX(i)  DZ(k) 


+  hj)ng-  DZ(k) 

♦  lUBDroW”  -  ^Ul)  Dy(  J)| 


Now,  one  more  pass  is  made  through  the  entire  grid,  this  time  omit¬ 
ting  the  momentum  equations  but  calculating  the  interface  velocities,  re¬ 
sulting  in  the  integration  of  the  internal  energy  to  time  (~) . 

X  -  T*1’  ♦  V "*»*)  DX(  1)  DZOO 

L  L  AMx£  L  2- 

+  Hy(  j)  Dz(k) 
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The  option  to  integrate  backwards  from  time  (~)  to  n  if  a  negative 
internal  energy  is  encountered,  is  not  available  in  this  version. 

2.2.4  PH2  Routine 

Here,  we  move  mass  across  the  fixed  boundaries.  Momenta  and  energy 
is  carried  across  with  this  mass  and  this  approximates  the  transport  terms 
omitted  from  the  momentum  and  energy  equations  in  PHI.  Please  refer  to 
Fig.  2  for  the  following  discussion. 

The  points  A,  B,  C,  D,  E,  F,  G,  and  H  are  the  eight  corners  of  the 
cell.  The  following  notation  will  be  followed:  side  AEFB  refers  to  the 
left,  side  ABCD  refers  to  the  top,  CDHG  refers  to  the  right,  BFCG  to  behind, 
EPGH  to  bottom,  and  ADEH  to  the  front. 

The  five  quantities  associated  with  each  interface  are  as  follows: 


TOP 


RIGHT 


BOTTOM 


LEFT 


BEHIND 

— 


AMPY  =  mass  crossing  the  top 

AMUT  =  X  momentum  of  this  mass 

AMVT  =  Y  momentum  of  this  mass 

AMWT  =  Z  momentum  of  this  mass 

DEIET  =  specific  energy  -across  the  top 

AMMP  =  mass  crossing  the  right 

AMUR  =  X  momentum  of  this  mass 

AMVR  =  Y  momentum  of  this  mass 

AMWR  =  Z  momentum  of  this  mass 

DEIER  =  specific  energy  across  the  right 

AMMY  =  mass  crossing  the  right 

AMMU  =  X  momentum  of  this  mass 

AMMV  =  Y  momentum  of  this  mass 

AMMW  =  Z  momentum  of  this  mass 

DELEB  =  specific  energy  across  the  bottom 

GAMC  (j)  =  mass  crossing  the  left 

FLEFT  (j)  =  X  momentum  of  this  mass 

YAMC  ( j)  =  Y  momentum  of  this  mass 

ZMOM  (j)  =  Z  momentum  of  this  mass 

SIGC  (j)  *  specific  energy  across  the  left 

BMASS  (M)  =  mass  crossing  the  back  surface 

BXMCM  (M)  =  X  momentum  of  this  mass 

BYMCM  (M)  =  Y  momentum  of  this  mass 

BZMCM  (M)  =  Z  momentum  of  this  mass 

BENR  (M)  =  specific  energy  across  the  back  surface 
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FMASS  «=  mass  crossing  the  front 
TO0NT  FXMOM  =  X  momentum  of  this  mass 

- - —  FYMCM  =  Y  momentum  of  this  mass 

FZMOM  =  Z  momentum  of  this  mass 

FENR  =  specific  energy  across  the  front 

Following  the  typical  cell  (L)  the  masses,  the  momentas,  and  the 
specific  energies  are  now  available  at  the  left  and  bottom  and  back  bound¬ 
aries  of  cell  (L)  from  the  previous  column  sweep,  the  cell  below  and  the 
previous  X-Y  slab  (K-l). 

The  proper  boundary  conditions  are  first  set  for  the  cell  (L®l).  We 
then  begin  by  calculating 


VABOVE 


V  +  V 

L  L+iMAX 


for  the  top 


\  +  \+1 

URR  =  - 5 -  for  the  right 


WOUT  for  in  front 


Then  form 


and 


VABOVE 


VABOVE 


1  +  <\+iMAX  '  V 


URR 


Ay(  J) 

URR 


At 


WOUT 


1  •  V  . . 

Ax(i) 

WOUT 


1  +  (V, 


f 

I/K  iMAX)(  JMAX) 
""  Az(kj 


-%) 


At 


Now  we  can  calculate  the  mass  crossing  the  three  boundaries  as 


AHpop  -  AMFY  = 

“right  -  "wp 


AMX(M)n  VABOVE  At 
”  Ay(  J) 

AMX(M)n  URR  At 
=  Ax(i) 
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_  AMX(M)  WOUT  At 

^RONT  =  mASS  =  -  ,V,k('k') - 

where  (M)  is  the  index  of  the  donor  cell.  The  donor  cell  is  calculated 
from  the  sign  of  the  weighted  velocity. 

The  momenta  associated  with  these  three  masses  are  now  calculated 
where  the  velociCj^  in  the  momenta  is  from  the  donor  cell.  The  total 
specific  energy  that  these  mass  fluxes  carry  are  calculated  at  this  time 
also.  The  momenta  associated  with  the  flux  at  the  top  is: 

X  component  =  AMUT  -  AMFY  [tf(N)] 

Y  component  =  AMVT  =  AMFY  [V(N)] 

Z  component  *  AMWT  =  AMFY  [ft(N)] 


Specific  energy  ■  DEIET  =  ^ - *■ 

Those  associated  with  the  flux  at  the  right  are: 

'i  U 

X  component  =  AMUR  =  AMMP  [tf(N)] 

Y  component  =  AMVR  =  AMMP  [V(  N)] 

Z  component  *  AMWR  =  AMMP  [V?(N)] 

[Tf  2  +  V  2 

Specific  energy  -  BEIER  =  7,,.,  .  SSl—JJ1) 


+  V(»)] 


Those  associated  with  the  flux  in  front  are: 

X  component  =  IMASS  [tf(N)]  =  FXMOM 

Y  component  =  K4ASS  [\T(U)]  ■  FYMOM 

Z  component  =  IMASS  [W(  N)]  ■  FZMOM 

Specific  energy  =  - - - *  • 

Where  again  (n)  «  index  of  the  donor  cell. 
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The  mass  now  in  cell  ( L)  is  equal  to  DELM  =  AMX(  L)  +  GAMC  ( j)  +  AMMY  - 
AMPY  -  AMMP  +  BMASS  M)  -  MASS  which  equals  the  original  mass  plus  the  mass 
flow  across  the  left,  the  bottom,  behind,  less  the  mass  flow  across  the  top, 
right  and  in  front. 

The  total  X  momenta  that  has  come  into  or  left  cell  ( L)  is  =  SIGMU  ' 

Fleft(  j)  +  AMMU  +BXMOM(M)  -  AMUT  -  AMUR  -  FXMQM  =  the  momenta  crossing  the 

left  boundary  plus  the  momenta  crossing  the  bottom  boundary  plus  the  momenta 
crossing  the  back  boundary  less  the  momenta  crossing  the  top,  the  right  and 
front  boundary. 

The  total  Y  momenta  that  has  come  into  or  left  cell  ( L)  =  SIGMV  = 

YAMC(j)  +  AMMV  +  BYMCM(M)  -  AMVT  -  AMVR  -  FYMOM  =  momenta  crossing  the 

left,  the  bottom  and  behind  less  the  momenta  crossing  the  top,  right  and 
front  boundary. 


The  total  Z  momenta  that  has  come  into  or  left  cell  (L)  is  *  SIGMW  * 
ZMOM( j)  +  AMMW  +  BZMQM(M)  -  AMWT  -  AMWR  -  FZMCM  =  momenta  crossing  the  left, 
the  bottom  and  behind  less  the  momenta  crossing  the  top,  right  and  front 
boundary. 

Now  by  conserving  momenta  and  total  energy,  we  can  calculate  new 
velocities  and  specific  internal  energy  of  cell  (L) 

«»!,  +  mue  +  MUBH  -  MUT  •  mur  -  Mlt  +  mu(l)  '  (DE“)  “(I) 

and 

“le  +  “b  +  ^BH  -  °T  -  +  m(L)  -  (DEIM)  va) 

and 

^LE  +  +  ^BH  ~  mT  -  mR  -  +  m(L)  =  ^  wn+l 

and  solve  for  the  three  velocities  at  cycle  n+1. 


The  new  specific  internal  energy  is  the 

.n+1  ELE  +  +  ET  "  ER  "  EF  +  EL 

(L)  = 


total  less  the  kinetic  = 

(u^1)8  -  <v"+1)2+  <wf  ^ 


DEIM 


2. 
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DELM  «=  the  new  mass.  The  subscripts  LE,  B,  BH,  T,  R,  F,  L,  refer  to 
the  left,  bottom,  behind,  top,  right,  front  and  cell  in  question. 

Now  the  five  variables  that  were  calculated  for  the  top  of  cell  (L) 
become  the  bottom  quantities  for  cell  (L  +  iMAX)  and  the  variables  that 
were  calculated  for  the  right  of  cell  (L)  become  the  left  quantities  for 
cell  (L+l)  and  the  variables  that  were  calculated  for  the  front  of  cell 
(L)  become  the  back  quantities  for  cell  L  +  ( iMAX ) ( jMAX ) .  This  completes 
the  calculation  for  cell  (L).  After  completion  of  PH2,  all  that  remains  to 
complete  cycle  n+1  is  to  update  the  pressures,  which  is  done  in  the  CDT 
routine. 

3.  TEST  PROBLEMS 

A  series  of  test  problems  were  undertaken  to  check  out  the  TRIOIL 

code. 

The  early  test  problem  consisted  of  8  cells  (2  in  the  x,  2  in  the 
y  and  2  in  the  z  direction)  containing  internal  energy  only,  free  to 
expand  into  a  vacuum.  The  purpose  was  to  check  the  free  surface  treatment, 
the  symmetry  of  disturbance  in  the  three  directions  and  the  possible  spher¬ 
icity  of  the  expansion  at  large  distances.  The  results  were  encouraging, 
exact  symmetry  existed  and  the  expansion  was  spherical  considering  the 
coarse  zoning. 

The  second  test  problem  (as  shown  in  Fig.  3)  consists  of  a  grid 
11  x  11  x  11.  The  corner  cell  has  lO1^  ergs/g,  surrounded  by  like  material 
but  cold.  The  three  adjacent  sides  were  treated  as  reflective  boundaries, 
while  the  three  opposite  sides  were  transmittive  boundaries.  A  check  of 
computer  results  were  made  against  the  G.  I.  Taylor  strong  shock-point 
source  solution.  The  comparison  is  presented  in  Fig.  4.  As  indicated, 
the  position  of  the  shock  front  is  in  good  agreement  with  theory,  however 
the  magnitude  is  somewhat  less  in  the  TRIOIL  solution.  The  latter  results 
are  consistent  with  similar  calculations  performed  with  the  OIL  (Ref.  1  )  code. 

The  third  test  problem  was  a  normal  impact  to  be  compared  with  the 
OIL  code.  The  configuration  (as  indicated  in  Fig.  5)  for  the  TRIOIL  code 
was  a  cube  of  4  x  4  x  4  cells  impacting  normal  to  a  semi- infinite  target. 


The  configuration  for  OIL  consisted  of  a  right  circular  cylinder  (D  =  L) 
where  the  diameter  was  equal  to  the  side  of  the  cube  in  TRIOIL.  The  OIL 
problem  had  2  radial  zones  and  4  axial  in  the  projectile.-  The  velocity  of 
the  projectile  was  2.6  x  10^  cm/sec,  and  the  projectile  and  target  were 
both  aluminum. 

The  momentum  and  times  of  the  two  problems  were  scaled  by  the  ratio 
of  their  masses,  the  distances  by  the  cube  root  of  the  ratio.  Figure  6 
presents  the  total  positive  y  momentum  (axial,  in  the  case  of  OIL)  versus 
time  for  the  2  codes.  The  agreement  is  excellent.  Figure  7  displays  the 
pressure  attenuation  into  the  target  (along  the  axis  for  OIL  and  in  either 
of  the  2  inner  original  columns  for  the  TRIOIL).  Here  again,  we  expect 
some  difference  at  early  times  since  one  projectile  is  a  cylinder  and  the 
other  a  cube,  but  the  agreement  is  very  good. 

Figure  8  is  a  velocity  (in  the  x-y  plane)  plot  for  the  plane  of  sym¬ 
metry,  z  =  0.  Figure  9  is  a  velocity  plot  (in  the  y-z  plane  at  a  value  of 
x  in  the  center  of  the  projectile.  Exact  symmetry  is  indicated  by  these 
two  plots.  Figure  10  is  at  a  later  time,  presenting  the  velocity  (ir.  the 
x-y  plane)  plot  for  the  plane  of  symmetry,  z  =  0. 

Figure  11  is  a  velocity  (in  the  x-y  plane)  plot  for  the  plane  of  sym¬ 
metry,  z  =  0.  for  an  oblique  impact  (45°).  The  dark  lines  Indicate  the 
position  of  the  original  projectile  to  target  at  time  t  =  0.  Figure  12  is 
a  velocity  plot  in  the  y-z  plane  at  a  value  of  x  at  the  right  hand  side 
of  the  original  projectile. 

Since  the  original  formulation  of  the  TRI0IL  code,  investigations 
into  the  effect  of  obliquity,  and  velocity,  for  semi-infinite  and  thin  targets 
have  been  completed  successfully. 
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Fig.  7 --Pressure  attenuation  into  the 
target  comparison  for  the  2  codes 
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Fig.  9--Velocities  (Y-Z)  in  the  Y-Z  plane  at  a  X  position 
of  the  original  center  of  projectile  for  90  degrees 
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Fig.  IO--Velociti.es  (X-Y)  in  the  plane 

of  symmetry  (Z=0)  for  90  degrees 
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Fig.  H--Velociti.es  (x-Y)  in  the  plane  of 
symmetry  (Z=0)  for  45  degrees 


OOOOOOOOOBOOO  oooooo 


FORTRAN  LISTING 


34 


4. 


1  WALLY 


S V  07-25-67  1720 


*****  note »  the  Following  set  of  dimension* 
common*  and  equivalence  is  to  bl  used  for  all 
subroutines  with  the  exception  of  subroutine  cards . 

DIMENSION  AIX(oOGO) f AMX (6000 ) fU(6000) fV(6000) *W(6000) fP(6000) f 
luX(3u) fUY(30) fD2(30) »UL(30) fPL(30) fX(30) fY(30) *ZCoR(30) f 
2Z( 150) f IZ ( 150 ) » FLLFT ( 30) » YAMC ( 30 ) »  SIGC (30) fGAMC ( 30) f2MOM(30) f 
3PdINU ( 700 )  *  UbIND  (700)  fBMASS(700)  *  BXMOM (700)  * BYMOM ( 700 )  f 
4i3ZM0M(70U)  fBENR(7Q0) 

DIMENSION  PR (50) fRK(30) 

COMMON  ZfAIXf AMX *  U  *  V » W f P , DX , DY  f DZ f 
IuLfPLfXf YfZCORf PK» SIfaC* GAMC» 

2ZM0M r  UbINQ  f DM ASS  *  BYMOM * 

3aZM0h fUENRf AREA* DIG * BOUNCE ' PABOVE f PBLOf 
4PIDTS  *  PRR  *  RHO  *  SI G*  UVMAX  *  V ABOVE  *  VBLO  * 

5VEL f  WPS f WS * WSA  t WSQ f WSC  *  I  *  1 1  *  IN  * IR  * 

6I.vS»lWSA»  IWSb*  IWSC* J»  JNfJPf  JR»K»KDT»KN» 

7KP » KR  >  KRM *  L  *  M »  M A  *  MB  * MC *  MD  *  ME  *  MZ * N 
COMMON  REZ  *  TRAD  *  DTRAD  f  RADEB  *  RADER  ,  RADET f  XI f  X2  f  Yl r  y2  r  IM.AXA 
EQUIVALENCE  (Z* IZ,PKOB) f (2(2) t CYCLE) f (2 (3) fDT) f 
1 (Z (4) *  PRINTS) * (Z(b) fPRlNTL) > (Z(6) » DUMPT7 ) * ( Z ( 7 ) fCsTOP) * (2(8) fPIDY) 
2* (Z(9) »6AM) * ( Z ( 10 ) fGAMD) * ( Z ( 1 1 ) t GAMX ) * (Z(12) fETH) , (2(13) fFFA) * 

3 ( Z ( 14 ) f  FFB ) * (Z(15) t TMASS) * (Z(l6) fXMAX)f (2(17) *  YMAX) * (Z(18) f2MAX) * 

4  (Z  ( 19 )  *  DNN )  t'  (2(20)  *DMIN)  *  (2(21)  *  DTNA )  *  (Z(22)  t  REZFCT )  *  (Z(23)  fTOZONE 

5)  , (Z(24) fECK) , (Z(25) fSBOUND) f (Z(?6) fCABLN) * (Z ( 27 ) , T ) f ( Z ( 28) *  GMAX ) f 
6 ( Z ( 29 )  * WSGD )  f  (Z(30)  fWSGX)  f  (Z(3l)  *GMADR)  *  (Z(32)  fGMaXR)  f  (Z(45)  fDTCHK 
7) * (Z(46) fPCSTAB) * (Z(47) fCNOT) f ( Z ( 46 ) f BFACT) * (Z(49) fEPSI) f (Z(50)»S1 

6)  j (Z(51) fS2) f (Z(52) »S3) f (Z(53) fS4) f (Z(54) fS5) r (Z(5S  fS6) » 

9(Z(5b) »S7) f (Z(57) fS6) # (Z(56) ;S9) f (Z(59) fSIO) 

EQUIVALENCE  (2(60) fAMLOST) f (Z(61) fELOST) f (Z(62) fXmLOST) f 
1(Z(63) fYMLOST) f ( Z  ( 64 ) fZMLOST) f (Z(65) *ENE6)r (Z(66) ,RHONOT) , 

2(2(6 7) fVELOC) f (Z (68) fBUG) f (Z ( Ql ) fNPR) f (Z(82) fNPRI ) f 
3(Z(63) *NC) * (Z(Q4) fNPC) * (Z(05) fNRC) f (Z(66) fIMAX) f (2 (87) f JMAX) f 
4(Z(88) fKMAX) f ( Z ( 89 ) » KM AX A) f (Z(9Q) f IXMAX) f (Z(91) fNOD) r 
5(Z(92) fNOPR) f (Z(93) fII) » (2(94) f 12) f (Z(95) fI3) f (Z(96) f 14) f 
6(2(97) »N1) * (4(98) fN2) * (2(99) *n3) f (Z(100) fN4>  f (2( 10l>  »N5) » 

7(Z( 102) fN6) * (2(103) *N7) f (2(104) f NO) r (2( 105) fN9) f (2(106) fNIO) f 
6(2(107) »N11) » (2(108) *K1) f (2 ( 109) »K2) * 

9(2(110) » Jl) f (2(111) fJ2) 

EQUIVALENCE  (BMASSfPBIND) . (QXMOMf UBIND) f 
KULfFLEFT)  f  (PLfYAmCfPK) 


FOR  MAIN/S  f  MAI  N/S*M.A  I N/SS 


INPU076Q 


*********  3DOIL  *************************** 


MAIN0050 

THE  INPUT  ROUTINE  WILL  READ  A  TRIOIL  BINARY  DUMP 

tape  or  will  call  for  subroutine  set-up 

WHICH  WILL  GENERATE  A  DATA  TAPE 
FOR  RESTRICTED  GEOMETRY 
CALL  INPUT 


MAIN006Q 


o  o  o  o  ooooo  oooono  oonoo  oon  o  r>  o  o  o  o  oo 
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V* 


•4 
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COT  ROUTINE  CALLS  FOR  THE  EQUATION  OF  STATE # 

CALCULATES  THE  DT (HYDRODYNAMICS  TIKE  STEP)  AND 
ADVANCES  THE  CYCLE  NUMBER 

10  CALL  CDT  MAIN0070 

IN  EDIT  #  DETERMINE  WHETHER  TO  EXECUTE  A 
SNORT  PRINT#  LONG  PRINT#  BINARY  TAPE  DUMP 
OR  TO  STOP  THE  PROBLEM 

call  edit  mainooqo 

***  SENSE  LITE  1  SIGNIFIES  THAT  THIS 
IS  THE  LAST  CYCLE  OF  THIS  RUN. 

THE  LITE  IS  TURNED  ON  IN  THE  EDIT  ROUTINE . 

CALL  SLITET(1#K0U0FX)  MAIN0090 

GO  TO(30»20) #K000FX  MAIN0100 

PHl  INTEGRATES  THE  MOMENTUM  AND 
ENERGY  EQUATIONS  DUE  TO 

TO  PRESSURE  FORCES  ONLY  #  THE  CONVECTIVE 
terms  are  Temporarily  dropped. 

20  CALL  PHI  MAINOHo 


PH2  SOLVES  ThE  MASS  TRANSPORT  EGUATIONS 
AND  MOVES  MOMENTUM  AND  ENERGY 
ACROSS  THE  FIXED  CELL  BOUNDARIES 
TO  APPROXIMATE  THE  CONVECTIVE  TERMS  THAT 
WERE  OMITTED  IN  PHl 

call  ph2 

ALL  CELL  QUANTITIES  »  EXCEPT  THE  PRESSURE 
HAVE  NOW  BEEN  ADVANCED  TO 
CYCLE  N+l 

GO  TO  10 
30  CALL  EXIT 
END 

FOR  cards/s»cards/s»cards/ss 
subroutine  CARDS 

DIMENSION  TABLE ( 1 ) » CARD (7) »LABLE ( 1 ) 

COMMON  TABLE 

A  2  IN  COLUMN  1#  ROUTINE  WILL  FIX  THE 
FLOATING  PT.  NO. 

A  1  IN  COLUMN  1#  MEANS  THIS  IS  LAST  CARD  TO 
READ  IN. 

EQUIVALENCE (TABLE (1) »LABLE(l) ) 

'WRITE  (6#  10) 

1  READ  (5  # 11 ) IEND » LOC  #  NUMWPC » (CARD ( I ) # 1  =  1 #  NUMWPC ) 
*RITE (b» 12) IEND#  LOC » NUMWPC  # (CARD(I) » 1=1# NUMWPC) 
DO  4  1=1# NUMWPC 

J=LOC+I-l 
I F ( I END-2 ) 2  #  5  #  2 
5  LA8LE(U)=IFIX(CARD(I) > 

GO  TO  4 

2  TABLE (U) =CARD ( I ) 

4  CONTINUE 

IF( I END-1) 1»3» 1 

3  RETURN 


MAIN0120 

MAIN0130 


MAIN0140 

MAIN0150 

MAIN0170 

CARD0010 

CARD0020 

CARD0030 


CARD0050 

CARD0080 

CARD0100 

CARD0110 

CARD0120 

CARD0130 

CARD0140 

CARD0150 

CARDOloO 

CARDC170 

CARD0180 


OOOOOOOiD  ’  *  OOOOO 
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CAR0019Q 


formats 

10  FORMAT (20H1  TR10IL  INPUT  CARDS///) 

11  FORMAT (11*15*11*  GP7E9.4 )  CARD0210 

12  FORMAT ( 1H  I4,I7*I3*1P7E14.6)  CARD0220 

END  CARD0230 

QZ  FOR  UNCLE/S  *  UNCLE/S *UNCLE/SS 

subroutine  uncle 

REWIND  N7 

SUBROUTINE  ((  UNCLE)))  IS  CALLED  WHENEVER 
A  COuhD  ERROR  IS  encountered 
IN  ANY  SUBROUTINE*  ITS  MAIN  FUNCTION 
IS  TO  PRINT  THE  CELL  QUANTITES  OUT  IN  THE  FORM 
OF  THE  NORMAL  LONG  PRINT 
NK=9U 

WRITE (6  *  8120 ) NC 

8120  FORMAT (1H0//6UH  AN  ERROR  HAS  OCCURRED  AND  SUBROUTINE  UNCLE  HAS  BEE 
IN  CALLED  AT  CYCLE15////) 

5000  .jRITe  ( 6  *  8116 )  PROS  *  NC *  T » DTNA  *  TRAD  *  DTRAD*NR*  N1 *  N2  * N3*  N4 
uO  1126  KK=K1*K2 
y.RITe  ( 6  *  9041 )  KK  *  ZCOR  (KK )  #DZ(KK) 

5004  UO  5U5G  1=11*12 
WS1=1. 

J=J2+1 

K=J2*IMAX+I+(KK-1)*IXMAX 

UO  5046  L=J1*J2 

J=J-1 

k=k-imax 

5012  IF(AMX(K) )5046*504o*5014 
5014  IF ( WS1 ) 5019  *  5019*  5016 
5016  WRITE (6* 8135) I*X(I)*DX(I) 

WSl=0o 

5019  v.S=Amx(K)/(DX(I)*DY(J)*DZ(KK)  ) 

»SC=p ( K ) *1 . E+4 

5018  rtKlTE(6*8l08) J»U(K) *V(K) »WSC*AMX(K) *WS*AIX(K) *W(K) *Y( J) 

5046  CONTINUE 
5050  CONTINUE 
1126  CONTINUE 
RETURN 

8108  FORMAT (I3*1X*1P8E12.5) 

81160FORMAT(6H1PROBLEM6X*5HCYClE9X*4HTIME13X*2HDT13X*4hTRAD11X*5HDTRAD1EDIT3380 
12X»2HNR6X»2HN14X»2HN24X,2HN34X*2HN4/(F7.1*Ill*2X*lP4El6.7*I10r2X*4EDlT3390 
216)  ) 

81350FORMAT<1H  ///4H  I  =  13* 6X * 6HX < I )  =F12 . 3 * 6X*7hDX ( I )  =F12.3//3H  J8X*EDIT3520 
11HX10X*1HY10X*  3HF/A9X*  3HAMX9X  *  3HRH08X  * 3HAIX9X  * 4H  W  8X*2H  Y/) 

9041  FORMAT ( 1H  /// 4H  K  =  13  * 6X * 9HZC0R ( K )  =F12.3*6X*7HDZ(K)  =F12.3) 

END 

I  FOR  sltup/s*setup/s*setup/ss 
subroutine  setup 

CALCULATE  THE  ADDITIONAL  INDICES  THAT  ARE  FUNCTIONS  OF 
I MAX  *  JMAX  r  AND  KMAX 

IXMAX= ( I MAX ) * ( JMAX) 

KMAXA=  KMAX*IXMAX 


SET  ALL  CELL  CENTERED  QUANTITIES  TO  ZERO 


onononnn  o  o  o  o  o  n  on  o  r>  oo  o  n 
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DO  1  K=1»KMAXA 
U(K)=Q. 

V(K)=0. 

WCK)=0. 

AIX(K)=0. 

AMX(K)=0. 

P(K)=0. 

1  continue 
X(1)=0X(1) 

CALCULATE  ALL  X,S 
DO  10  1—2 , 1 MAX 
X(I)=X(I-1)+DX(1) 

NOTE,  dx  is  CONSTANT  for  all  I 
0X(I)=DX<1) 

10  CONTINUE 
Y(1)=UYC1) 

CALCULATE  ALL  Y»S 
DO  11  J=2,JMAX 
Y(J)=YCJ-1)+DY(1) 

NOTE,  DY  IS  CONSTANT  FOR  ALL  J 
DY(J)=DY(1) 

11  continue 
2C0RC1)=DZ(1) 

CALCULATE  ALL  ZCOR,S 

DO  12  K=2,KMAX 

ZCOR  C  K ) =ZCOR ( K-l ) +DZ ( 1 ) 

NOTE»  DZ  IS  CONSTANT  FOR  ALL  K 
DZ(K)=DZ(1) 

12  CONTINUE 

rhonot  is  initial  desity  for  all  material 

J3=S1 

DO  100  K=1»KMAX 
LL=(K-1)*IXMAX 
DO  100  I=1»IMAX 
DO  100  J=J3,JMAX 
L=LL+(J-1)*IMAX+I 
AMX (L )=DX( I )*DY(J)*DZ(K) ♦RHONOT 
100  CONTINUE 

si=interface(j)  value  +1  between  projectile  and  target 

S2=BACK  BOUNDARY  +1  OF  THE  PROJECTILE ( K )  VALUE 

S3=  FRONT  BOUNDARY  OF  PROJECTILE (K >  VALUE 

S4=  LEFT  BOUNDARY(I)  VALUE  +1  OF  THE  PROJECTILE 

S5=  RIGHT  BOUNDARY  (I)  OF  THE  PROJECTILE 

S6=  BOTTOM  BOUNDARY  (J)  +1  OF  THE  PROJECTILE 

S7=  TOP  BOUNDARY (J)  OF  THE  PROJECTILE 

ETH=0. 

K11=S2 

K22=S3 

I11=S4 


.  DO  200  K=Kll»K22 

LL=(K-1)*IXMAX 
DO  200  1 = I 1 1 f 122 
00  200  J=Jll»J22 
L=LL+(J-l)*lMAX+I 
AMX(L)=DX(I)*DYCJ)*DZ(K)*RHONOT 
C  VEL0C=1NITIAL  Y  COMPONENT  OF  VELOCITY 
V (L)=VLLOC 

C  S9=  INITIAL  X  COMPONENT  OF  VELOCITY 

U(L)=S9 

C  S10=  INITIAL  Z  COMPONENT  OF  VELOCITY 

w(L)=SlO 

ETH=ETH+AMX(L)*(U(L) **2+V(L) **2+W(L)**2)/2. 

200  continue 
C 

C  PRINT  THE  QUANTITIES  ASSOCIATED  WITH  THE  GRID 
C 

WRITE (6 , 8000) IMAX, ( X ( I ) » 1  =  1# IMAX) 
*RI7E(6»8003)IMAX# (DX(I) »I=1»IMAX) 

WRITE ( 6 » 8001 ) JMAX » (Y(J) # J=1»JMAX) 

WR ITl ( 6 » 8004 ) JMAX , (DY(J)  #J=1»JMAX) 
wRl TE (6  f  8002) KMAX # (2C0R(K) #K=1#KMAX) 

WRITE ( 6#  8005) KMAX  # (DZ(K) ,K=1»KMAX) 

WRITE (6* 8006) IMAX , JMAX  » KMAX » IXMAX#KMAXA 

8000  FORMAT ( 1H  /10H  X ( I ^  1=1 » 12/ ( 5F16.6) ) 

8001  format (in  /ioh  y(j)  j=i# i2/<5Fi6.6> > 

8002  FORMAT ( 1H  /13H  ZCcR(K)  K=1 *  12/ (5F16.6) ) 

8003  FORMATUH  /11H  DX  { I )  1=1 » 12/ (5F16.6) ) 

8004  FORMAT ( 1H  /11H  DY{J)  J=1 » I2/(5F16.6) ) 

8005  FORMAT ( 1H  /11H  DZ(K)  K=1 , 12/ (5F16.6) ) 

8006  FORMAT (718) 

C 

C  k.RlTE  A  DUMP  TAPE  (FOR  T=0.)  FOR 
C  THE  TRIOIL  CODE 

REWIND  N7 
WS=555.0 

wRl TE ( N7) WS » CYCLE » N3 
WRITE ( N7 ) (Z(I) »I=1#150) 

.vRITe  ( N7)  (U  ( I )  » V  ( I )  #  W  ( I )  , AMX  ( I ) » AIX ( I ) » 1=1 # KMAX A) 
WRITE (N7) <X(I) » 1=1# IMAX) 

WRITE (N7) (Y(J) , J=i# JMAX) 

WRITE (N7) (ZCOR(K) ,K=1#KMAX) 

WS— 666 • 

WRITE(N7)WS»WS»WS 

return 


end 

(01  FOR  INPUT/S# INPUT/S# INPUT/SS 


subroutine  input 

iNPUOOlo 

c 

MAIN0020 

c 

INPU0900 

CALL  SLITE  (3) 

INPU0980 

c 

IWSA=THE  NUMBER  of  BCD  HEADER  CARDS  TO 

c 

READ  IN 

READ (5*8007) IWSA 

noon  . .  noon  rvn  n 
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□0  7  I=1*IWSA 

HEAD  (b»8004)lWS  INPU1000 

WRITE (6*6004)  IWS 
.  *  7  CONTINUE 

•  6  CALL  CARDS  INPU1020 

NOTE*  PROVSIONS  FOR  CALLING  A  GENERATOR 

CODE  SUCH  AS  (CLAM)  »  DOES  NOT  EXIST  IN 
3  DIMENSIONS  AT  THIS  TIME 
IF (PK ( 3) ) 8687 * 8888*8888 
8888  CALL  CARDS 
CALL  SETUP 


8887  CONTINUE  INPUIO60 

GO  TO  1000  INPU1090 

10  CONTINUE  INPU1120 

CALL  CARDS  INPU1130 

GO  To  2000  INPU1140 

40  uO  4b  K=1»KMAXA  INPU1160 

45  P(K)=0.0  INPU1170 

integrate (the  time  and  cycle  number  )  backwards*  since  they 

WILL  BE  ADVANCED  IN  CDT  SUBROUTINE 

T=T-uTNA  INPU1180 

NC=NC-1  INPU1190 

CYCLE=NC  INPU1200 

NPC=NPC-1  INPU1210 

UVMAX-0  •  0  INPU1220 


CALCULATE  DX*DY*DZ  SINCE  THEY  ARE  NOT  ON  THE  TAPE 


’aS=0, 

DO  50  1  =  1 » IMAX  INPU1230 

DX(I)=X(I)-WS 

WS=X(I) 

50  CONTINUE 

ftS=0. 

DO  5b  J=1*JMAX  INPU1250 

UY(J)=Y(J)-WS 

WS=Y(J) 

55  CONTINUE 
WS=0. 

DO  65  K=1*KMAX 
UZ (K) rZCOR (K) -WS 
WS=ZCOR(K) 

65  CONTINUE 

DUMP  THE  Z  BLOCK 
K=1 

DO  80  1=1*9 
L=K+b 

.'•'RITE (6* 8005) K »  (Z(N)  *N=K»L) 

80  K=L+1 
K=61 

DO  81  1=1*8 
L=K+a 

WRITl(6»8006)K* (IZ(N) »N=K»L) 

81  K=L+1 
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GO  To  10000 
C 
C 

.C  READ  THE  til  NARY  TAPE  FOR  A  RESTART 
/•* 

V- 

1000  MZ=loQ 
i  *-.;S=u 

lu03  aE'aIU'J  N7 

100 a  REmD(u7)Pk(1)  *PR(2>  *N3 

i\!R=NG  +  o 

100b  IK (PR (15-550.0) 1010 » 1016* 1010 
lulO  1».S-1'aS+1 

lull  if  IM0J<  1aS*3)  >9902*9902*1003 
101b  IF  (P|<(2))l0l0*101c*10ia 
1016  Ir  U-'K(2)-PR(2> ) 1023* 1023*1020 
C  RcAD  OVER  *  Tj-,1  S  IS  f JOT  THE  CORRECT  CYCLE 
1020  GO  1u2E  L=2*NR 

1022  READ  (  n'/  )  V<bS 
oO  Tu  1 0  U  4 

1023  UEAD(u7) (Z(l) *  I  =  1*M?J 

c  check  for  correct  problem  numuer 

IF ( AbS (PRC5-PK ( 1 ) 01) 1024,1024*9901 

1024  KtAD(N7> ( U ( I ) ?  V ( 1 ) *  W ( I ) *  AMX ( I ) *  A I X  C I ) , I=1*KMAXA) 
REAL)  ( u7 )  (  x  C I )  *  1  =  1.  I  MAX  ) 

READ  ( r7 )  (Y(J)  * J=i* JMAX) 

I<EA0(r7)  (2COK(K>  *  K-l » KMAX ) 

1034  READ ( r7 ) PR \ i ) ,PK(2) *PK(3) 

1036  IF ( PR ( 1 ) -bb5. 0 ) 9904  *1040  *  1036 
1036  lF(PR(2)-666r0)99o5*lG40*9905 
1040  oO  To  10 

2000  IF('a,EGX)9906,20iG,2005 
C  CALCULATE  l./( GAMMA-1.  )  AND  GAMMA/ ( GAMMA-1 . ) 
2uGb  GAMX=i . 0/ ( WSGX-1 . 0 ) 

2010  i\'6GX=  ( GAMX+i  •  0  )  / GAMa 
G  M  A  X  R — G  A  M  X  *  WSGX 
2012  IF  C.VbGD)  9907, 2020, 2015 
2015  GAMD=1.0/(w5bD-l«0) 

2020  ..'SGD=  vCAMD+i .  0  )  /GAME 
vjMADR=GAMLMWSGQ 
GMAX=w6GD 

C  SEARCH  FOR  MAX  GAMMA 

IF (WSGD-wSGX) 2025, 2030 *2030 
2025  GMAX=WSGX 
2030  GO  To  40 
C  ERROR 

9901  hk=1023 

GO  TO  9999 

9902  HK— lOll 

GO  To  9999 

9904  uK=l03o 

GO  To  9999 

9905  N*=1G36 

bO  To  9999 

9906  NK=2000 

uO  TO  9999 

9907  NK=20l2 
9999  NK=1 


INPU1370 

INPU13B0 

INPU1390 


INPU1400 

IMPU1410 

INPU1420 


INKU1460 
INPU147  0 

INPU1460 

INPU1490 

INPU1500 

INPU1510 

INPU1530 


INPU1550 


INPU1680 

INPU1690 

INPU1700 

INPU1750 

INPU1760 

INPU1770 

IKPU1780 

INPU1790 

INPU1800 

INPU1810 

INPU1820 

INPU1830 

INPU1040 

INPU1850 

INPUI80O 

INPU1900 

INPU1910 

INPU1920 

INPU1930 

INPU1940 

INPU1950 

INPU1960 

INPU1970 

INPU19S0 

lNPU19q0 

INPU2000 

INPU2010 

INPU2020 


oooooooo  ooo  oooo.  ooo 


WRITE  ( t>»30 00)  PR(1)  ,PR(2)  , PK ( 1 ) »PK (2) *PK (3) # IWS 
CALL  UNCLE 
CALL  DUMP 


4l 


C.  •  INPU2040 

10000  RETURN  INPU2050 

C  INPU2060 

C  •  FORMATS  INPU2070 

8000  FORMAT ( 1P5E14.6, 15) 

80040FORMATCIl»7lh  INPU2090 

1  )  INPU2100 

300b  FORMAT ( 14, 1a» 1P9E12.5) 

3006  FORMATC 14, 1X,9I7) 

3007  FORMAT  C613) 

C  INPU2120 

END  1NPU2130 

CZ  FOR  CUT/S, CDT/S*CDT/SS 

subroutine  cut  cdt  0010 

c  CDT  0020 

3000  VLL=0.0  CDT  1030 


set  up  the  loops  for  calculating  The  material  pressure 

DO  3050  K=K1»K2 
LL=(K-1)*IXMAX 
3005  DO  3050  1=11,12 
3010  L=I+LL+(Jl-l)*IMAx 
3015  00  3050  J=J1»J2 
3020  lF(AMXCL) )9901, 3050, 3025 

THE  ES  ROUTINE  CALCULATES  THE  PRESSURE,  FOR  INPUT  IT  NEEDS  I,J,K,L 
AND  THE  AIX (ENERGY)  AND  AMX(MASS) 


3025  CALL  ES  CDT  1110 

THE  OUTPUT  FROM  ES  IS  THE  P( PRESSURE)  AND  MAX  (GAMMA-1.) 

3030  IF ( AbS l P (L ) ) -1 .E“20) 3035, 3035 , 3040 
3035  P ( L ) =o • 

3040  IF (WSGX”VEL)3050, 3050, 3045  CDT  1150 

5046  VLL=HSoX  CDT  1160 

3050  L=L+1MAX 

3055  KUT=1  CDT  1180 

U\/MAX=-1.0  CDT  1190 


now  set  up  the  loop  for  calculating  dt  from  the  particle 
velocities  and  the  courant  condition 

WE  FLAG  the  cell  that  is  controlling  the  time  step 
and  STORE  THE  VALUES  OF  I,J»K  INTO 
NIG ' Nil  *  AND  N9  RESPECTIVELY. 

DO  325b  K=K1,K2 
t-L=  (K-l )  *  IXMAX 
3070  DO  3255  1=11,12 
3075  L=I+LL+(U1-1)*IMAX 
3095  DO  3255  U=J1,J2 
3100  CONTINUE 


o  o  o  o  o  o  o  o 
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j.r!  A|..x(L)  }49ul »  3265*4 

c 

IS  i'Hi'j  16  LE66  TiJ.-t'  ol'CH!'.  OMIT  T||:I  3T ABILITY  CjIcCK  FOR  THIS 

CuLL 

4 

IF  ( .a;.,X  <  u)  /  ( JX  {  i  >  N/.'  (  J  J  *02  ( K )  ;  -DT CHK  )  3266 1  32 66 » 3115 

31 1  S 

S  X  (  i  ) 

CLT 

126o 

3120 

»r  ( o Y (J ) —SlG ) 3i2s » 3130  *  3 13 J 

CDT 

1270 

312o 

S  i  0  —  L<  Y  v  J ) 

CDT 

1260 

3-30 

IP  '  0*.  ( K)  -SIC )  s!31 1  6130 » bistf 

•  61  SI 

w-  1  N;  *■  W  ,£  \  l\  ) 

C 

i-LAl  f-  CALCULATE  The  SPEED  OF  SOUND  »  THE  PERFECT  6AS  SPEED 

0^ 

C 

jw.  „  0  \  D,3/'DR#iO  FKC'-i  THE  MlT-L  ZCUATION  OF  STATE 

6 1  S  o 

-  r*  vC..UT)40n  J  r  i  0  J  o  #  *r  0  0 1 

4000 

..0-S..K',  (G;«  AX»r  (L)/Aiv:X(L)  f  DX  (I )  >DY  £  U )  *i)Z  £ K  )  ) 

uO  o  32  Ob 

CDT 

1310 

4uOl 

i<  6  «  —  A '  S  ( p  (  L  )  ;  *  I  » !_  *r  4 

»'-,4Cj;oT  +  jF/\C  t  *  ;usa**EPSI) 

..  —  VS'!'  -1  •  L.  ”  W 

CUT 

134Q 

3i.(;6 

s  -  w  s  /  S  i  o 

CDT 

1360 

3210 

a  F  D,i  v  M  4  X  “  iii  6  )  Ocls*>j220*  j  e  2  U 

CDT 

1360 

321b 

I  'l  i  l 1  ~ 

CUT 

1370 

N  ii‘E 

—  i.' 

CDT 

1360 

u  y  —  p. 

UONA/.-.'.S 

CDT 

1390 

3„20 

oont.nu^ 

• 

W6-At/J{U(L)  >  /DX  ( 1 } 

322s 

IF  ( UVMA/>toS ) 3230  ' 323 s  *  323b 

CDT 

1460 

'  3^30 

U  V  iv'i  A  >;  Z  <h  $ 

CDT 

14  no 

rJl  0”. 

CDT 

1470 

N\)  =  K 

* 

Nl  1  —  U 

CDT 

1480 

3<-;3S 

io~  An.-.  i\l  ( L )  3  /  u Y  ( J ) 

3 1.4  0 

IF  (UVF.AX-wS )  324  6  e  3260  >  3260 

CDT 

1500 

Oc*+b 

,  ^  i  3  ~  i 

CDT 

1510 

i-j1  i  =  J 

CDT 

1520 

!J  J=K 

sir  if  a  —  i » b 

CDT 

1530 

3260 

CON  i  liROl 

CDT 

1540 

'«  u  —  A  Lj  6  ^  W  ( L  )  '  /  j  2  ( K  ) 

Ir  <  UVMAX-'jS)  6246  »  6250  »  5250 

6246 

N9  =  K 

NiO  =  i 
Nil-  j 
UVMAx-aS 
5260  o-ONTInUT. 

3266  L-l_+iM/'X 

IK  (UVilAX)  9912?  99X2  r  32o0 

here, check  the  o  pgssiule  options  fop  calculating  thf 

NL«;  iJ  T 

3260  ir  » C/,jLN) 90 r  91  * 3300  CDT  1560 

HcUL  OPTIONS  La  1ST  FOP  STARTING  THL  PHOuLLi' 
him  SMALL  TIPiL  STEPS  (A  SMALL  FRACTION  Op  STABILITY) 

FOi<  PROfjLLMS  PhLf'L  THL  INITIAL  ENERGY  IS 
r’NiMAKli-'r  i N  *  w.  Ki *lrw-  •  *  ■>  v  «  •<  •  *  «*  *  *  • 

90  IF  ( 2  ( o9}  -1 . 0  >  0  j  >  9«t  /  04 
90  2  ( 69 )  =2 (69 > -till  70 
60  TO  9b 


ooooc'i-onoooooo 


^3 

94 

2(b9)=1.0 

95 

UT=.o/VEL/UVMAX*PCSTAG*Z169) 

00  To  3295 

CDT 

1580 

•  •  91 

wS=UVMAX*  jT 

CDT 

159c 

• 

»sa=g.s/vll 

CDT 

1600 

3<;o5 

IF (FFA-WSA) 3276 #3276# 3270 

CDT 

1610 

.3270 

FFA  =  r.  bA 

CDT 

lt>20 

•3276 

IF ( W5-PFA) 3205#  33u0  #  3280 

CDT 

1630 

32oU 

uT=DT/WS*FFd/0.9 

CDT 

1640 

uO  Tu  3295 

CDT 

1650 

3265 

I F ( WS-P Fd ) 3290  #  3290  #  3300 

CDT 

1660 

3290 

UT=DT  *FFA/WS*0 .9 

CDT 

167c 

3295 

KDT=0 

CDT 

1680 

3300 

T=T+oTNA 

CDT 

1690 

62 

NC=NC+1 

CDT 

1780 

cycll=nc 

uvmax=uvmax*dt 

CDT 

1790 

NPC-NPC+1 

CDT 

1800 

3305 

IF  ( T) 9909#  3320  #3310 

CDT 

1810 

3310 

ir (KuT)9910#3315#3320 

CDT 

1820 

3315 

WRITE (6# 8000) T»DTnA#DT 

3320 

dtna=dt 

CDT 

1840 

GO  TC  3325 

CDT 

1850 

9901 

NK=3020 

CDT 

1870 

O0  To  9999 

CDT 

1860 

9909 

NK=350S 

CDT 

1890 

GO  TO  9999 

CDT 

1900 

.‘9910 

NK= 33 10 

CDT 

1910 

GO  TO  9999 

CDT 

1920 

9912 

NK-1 

GO  TO  9999 

9999 

NR=2 

WK I TE 1 6 »  60 02 ) I » J » K  #  L »  NK  #  NR »  NC 
WRlTE(o»8001) UVMAX  t AMX ( L ) »P(L) »DT »  VEL 

CDT 

1940 

6001 

FORMAT (1P5E14. 6) 

8002 

FORMAT (81 5) 

CALL  UNCLE 

CALL  DUMP 

3325 

return 

CDT 

I960 

8000 

FORMAT  (17H0CHANGE  DT  ...  T=1PE10.3# IcH 

DT(N)=1PE10.3#12H 

DT 

1 ( N+l ) =1PL1Q . 3) 

END 

CDT 

1990 

QX  FOR  PHl/S»PHl/S»PHl/SS 

SUBROUTINE  Phi  PHI  0010 

THE  VELOCITIES# ENERGIES#  PRESSURE  AMD  MASS  ARE  AT 
CELL  CENTERS 

FIRST  PASS  » CALCULATE  U,V  AND  W  FOR  CYCLE  N+l  #  AND  THE  WORK 
TERMS  USING  THE  VELOCITIES  AT  CYCLE  (N) . 

SECOND  PASS#  CALCULATE  THE  CONTRIBUTION  TO  THE  CHANGE  IN 
INTERNAL  energy  from  work  terms  evaluated  from  u#v  and  w 

AT  CYCLE  N+l 
2  PASSES  REQUIRED 


remember#  we  are  not  advancing  the  velocities  and  energy  to  CYCLE 

N+l# SINCE  WE  HAVE  NOT  AS  YET  SOLVED  THE  TRANSPORT  TERMS#  AS 
USUAL#  WE  REFER  To  THE  VELOCITIES  AND  ENERGY  FROM  THIS  ROUTINE 


C  Ab  T^LDA  QUANTITIES 

c 

C  CtLL  IN  QUESTION  =  L= ( J-l ) *IKAX+I+  (K-l) *IMAX*JMAX 
C 

C  IXMAa=IMAX*UMAX 

C 

C  KMAXA= ( IMAX) ( JMAXJ CKMAX) 

C 

C  CcLL  TO  THE  HlGHTz  L+l 

C 

C  CELL  ABOVE  =N=  L+IKAX 

C 

C  CELL  IN  FRONT  =NN  =L+IXMAX 

C 

c  n:  =flag  at  the  left 

C  N2  =F I  AG  AT  THE  RIGHT 

C  N3  =FLAG  at  THE  TgP 

C  l\4  =FLAG  AT  TrlE  BOTTOM 

C  N5  =FLAG  AT  THE  BEHIND 

C  i^o  =FLAG  AT  THE  IN  FRONT 

C 

C  SET  the  flags  for  INCREASING  OR  DECREASING  active  grid 
C  Counters  TO  ZEhO 

1X1  =  0 
1X2=0 
U  Y  1  =  0 
JY2=0 
KZ1=G 
KZ2=0 

C  SET  FLAG  FOR  SUBCYCLING 
VcL=i , 

C 

c  return  here  for  the  second  pass 

c 

C  SET  UP  THE  K  LOOP  FIRST 
2  DO  3  K=K1 t K2 
7  L= ( K-l ) * ( IXMAX ) +1 
C 

C  CHECK  FOR  BOUNDARY  CONDITIONS  AT  THE  LEFT 
C 

o  DO  330 2  J=1»JMAX 
6  IF (Ni ) 9  *  99 » 9 
C 
C 

c  tkansmxttive  left  boundary 

c 

99  PL ( J) =0 • 

GL(J)=U(L) 

GO  TO  10 
C 

c  reflective  left  boundary 

c 

9  PL(J)=P(L) 

UE( J) =0. 

10  L=L+I,v,AX 
3302  CONTINUE 

11  IF (K-l) 8999 » 12 » 7001 


o  o  ooo  ooo  ooo  ooo  ooo  o  r>  o  o  r>  r>  ooo  o  o  o  o 


^5 


7001  1F(K-KU15#12»15 

BACK  SIDE  QC.  HAVE  ALREADY  BEEN  SET 


12  CONTINUE 

BACK  SIDE  IS  TRANS#  OUT  WILL  TAKE  CARE  OF  IT  LATER 
uACK  SIDE  IS  REFLECTIVE 

13  UO  2302  N=1 > IXMAX 
IF(N5)6010»6011»6010 

6011  Pt3lND(N)=0. 

Ut3lNU(N)=W(N) 

GO  TO  2302 
6010  PdlNU(N)=P(N) 

UBIND(N)=0« 

2302  CONTINUE 

15  LL=(K-1)*IXMAX 

SET  DO  LOOP  IN  X  DIRECTION 

DO  4  1=11 # 12 

16  M=( J1-1)*IMAX+I 

17  L=LL+I+(J1-1)*IMAX 

SET  DO  LOOP  IN  Y  DIRECTION 

18  UO  5  J=J1»J2 
NN=L+IXMAX 

19  N=L+IMAX 

20  IF' J- 1)99 02? 21*7003 
7003  IF(Jl-vJ)3305#23»3305 

WE  HAVE  already  CALCULATED  THE  BOTTOM  BC. 
check,  for  bottom  boundary 

21  IF ( N4 ) 22 » 23  #22 

BOTTOM  BOUNDARY  IS  TRANS 
23  V8L0=V(L) 

pblo=o. 

GO  TO  3305 

BOTTOM  BOUNDARY  IS  REFLECTIVE 

22  VBLO=0. 

PbLO=P(L) 

NOW » 'WE  HAVE  THE  LEFT  BC,  (IF  REFLECTIVE)  AND  THE  BOTTOM ( IF  REFLECT) 

3305  IF(AMX<L) ) 9900 » 3340 » 3306 

CELL  IN  QUESTION  IS  VOID#  GET  OUT  AND  CONTINUE  THE  LOOP 

3306  IF ( IMAX-I ) 9901 » 331 1 » 3307 

WE  ARE  AT  THE  RIGHT  BOUNDARY  OF  THE  GRID 


46 

C 

C  Till-  TOP  IS  KU'LEClIVIi 

C 

30  PAUOVL=P(L) 

V  AI)OVl  =  0. 

00  TO  332M 

C 

C  THL  TOP  IS  r RAMS 

r 

31  PADOVL=POLO 
C 

c  Mou if- y  tin  For  Tuan*,  boundary 
c 

ETH=LTH-PAU0VE/2.4V  (L)  *DTtI)X(  I  >  *DZ<K> 

00  To  3323 
C 

C  WE  Akl  mot  at  tul  top 
3320  IK ( Al-iX  (N)  )  990b » 332? » 3324 
C 

C  CLLL  AIjOVL  IS  VOID 

3322  PAHOVL- 0 • 

3323  VAUOVL=V(L) 
oO  To  3326 

C 

C  NORMAL  FLOW  FOR  ALL.  CELLS  OCCUPIED 

3324  PAUOVL- ( P  ( L )  4  P  ( N ) )/2. 

32  IK  (  1  -  J )  332b »  33 »  99ii(> 

C 

C  DOTTOM  ROUNDARY  HAS  UCEN  SET 

C  •  WE  Are  AT  TUL  DOHOM 

33  IK  (N4  )  332b »  7UIJ0  >  332b 
C 

C  REFLECT  I VL  DOTTOM  bOUIJHARY  CONDITION  HAS  ALREADY 

C  ULLIJ  SET. 

C 

7000  PULO-pAbOVL 

LTH=l.Tli+PbL0/2.*V<L)'»DT*UX<I  )*DZ(K> 

332b  VADOVE=<V(L)+V(N) )/2. 

C 

C  CHECK  THL  l  DIRECTION 

C 

332M  IF (KMAX-K)99o7f 44lbf 4420 
C 

c  wl  all  at  tiil  front  ok  the  orid 

C 

44  ID  IK  ( I  Jo )  2999 1 34  1 2999 
C 

C  FRONT  lb  REFLECTIVE 

2999  PZH-P ( L ) 

WZR=0 . 

00  To  432h 
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C 

C  FRONT  IS  TRANS 

34  PZK=P|UNU(M> 

C 

C  MOUIf-Y  LTH  KOR  TrANS  BOUNDARY 

4419  LTII=LTH-PZK/2.  *  W  <L)  ♦l)T  H>X  ( I ) *DY ( J) 

00  TO  4323 

C 

C  CHECK  CELL  IN  FRONT 

C 

4420  IF‘(A|w,X(NN)  )  9900, 3;-,, 4324 
C 

C  CELL  IN  FRONT  '  lb  OCCUPIED 

C  CELL  IN  FRONT  IS  VOID 

C 

35  PZK=0. 

4323  WZK=W(L> 

00  TO  4320 

C 

c  normal  flow  in  2  DIRECTION  for  all  cllls  occupied 

4324  P2K=(P(LHP(NN)  )/2. 

IF  ( 1”K  )  432b »  37  » 99 (j 9 

C 

C  DC.  ULMINU  HAVE  ALREADY  BEEN  SET 

C 

C  WE  Are  IN  the  FIRST  (X-Y  ) PLANE  K=1 
37  lF(Nb)432b»0000,432b 
C 

C  REFLECTIVE  uc.  BEHIND  HAVE  ALREADY  BEEN  SET 
C  TRANS  UC.  IN  THE  BACK 
BOOO  PUINU{M)=PZK 
C 

C  MODIFY  cth  for  trams  boundary 

LTH=LTE'PZR/2.*W(L)*U.T+0Y(  J)  +DX  ( I ) 

432b  W2R=(W(L)+W(NN) )/2. 

c 

c  check  for  first  or  second  pass 
c 

4320  IF(Vt_L)9910»42»34u0 

c 

C  THIS  is  the  second  pass, skip  the  momenta  EQUATIONS 
r 
c 
c 

c  integrate  the  y  component  of  velocity  (V) 

3400  V(L)=V(L) +(PuEO-PABOVC)/AMX(L)*DT*DX( I)*DZ(K) 


on  oooo  on  oo  r>  n  n  ooo  ooo  oor>o  oooo 
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c  check  the  mass  to  the  right 

3307  IK(AMX(L+1> ) 9003 * 3312  *  33 14 

this  is  the  lsc.  at  the  right  of  a  occupied  cell  *  with  the 

NEIGHBOR  void. 

3312  PRR=U. 

3313  URR=U(L) 

GO  To  3310 

HERE*  VvL  ARE  AT  THE  RIGHT  BOUNDARY  OF  GRID  CI  =  lMAX) 

check  here  For  reflect  or  traijs  bc. 

3311  IF(M2)33ud*3309*3308 

REFLECTIVE 

3300  PRR=f*(L) 

URR=0. 

00  To  331o 

TRANSMITT  I VE 

3309  PRR=PL(J) 

MODIFY  LTH  HERE  AT  THE  TRAMS  BOUNDARY 

ETH=LTH-PRR/2.+U(L)*DT«'DY(  J)*DZ(K) 

GO  TO  3313 

HERE  IS  NORMAL  FLoW  FOR  ALL  CELLS  OCCUPIED 

3314  PRR=(P(L)+P(L*1) )/2. 

UKR=(U(L)+U(L+1) )/2. 

CHECK  HERL  FOR  ALONG  THE  LEFT  BOUNDARY  (1=1)  FOR  TRANS 
3316  IF ( I -1 ) 9911  *  SO  *  331 0 

80  IF(Nl)3310*bl, 331(i 

REFLECTIVE  (BUT  UC.  HAVE  ALREADY  BEEN  SET) 

TRANSM1TT1VE 

81  PL ( J) =FRR 

MODIFY  ETI I  AT  TRArS  BOUNDARY 

82  tTH=LlH+PRR/2.*U(L)*UT*DY( J)*DZ(K) 

00  TO  3310 

3310  IF (JMaX-J) 9904*3310 *3320 
C 

C  WE  ARE  AT  THL  TOP  OF  THE  GRID 
3310  IF ( N3) 30 » 31  * 30 


ooo  non  oooo  on’..  o.o 


IF(ABS(V(L) )-l.E-b) 3401 #3401, 3402 

3401  V(L)=0. 

INTEGRATE  THE  X  COMPONENT  OF  VELOCITY  (U) 

3402  U(L)=U(L)  +  (PL( J)-PRR)/AMX(L)*DT*DY( J) *DZ  UO 

40  IF ( AbS (U ( L) )-l.E-8) 3403 r3403» 3404 

3403  U(L)=0. 

integrate  the  z  component  of  velocity  (w) 

3404  h'(L)=W(L)  +  (PBIND(M)-PZR)/AMX(L)*DT*DY(J)*DX(I) 

41  IF(ABS(W(L) )-l.E-b) 3405,3405# 42 

3405  W(L)=0. 


HERE  CALCULATE  THE  CHANGE  IN  INTERNAL  ENERGY  DUE  TO  THE 
WORK  TERMS 

42  WS=P ( L) +UT/AMX (L) * ( ( VSLO-VABOVE ) /2 . *DX ( I ) *DZ(K ) 
1+(UL(J)-URR)/2.*DY(J)*QZ(K) 
2+(UBINU(M)-WZR)/2.*DX(I)*DY(J) ) 

43  AIX(L)=AIX(L)+WS 

CHECK  ON  ADVANCING  OR  DECREASING  GRID  COUNTERS 

5600  IF ( I-I2) 5999  ,5801,5801 

5801  IF ( 1X2) 5999  #5802#  5999 

5802  IF(ABS(UIL)  )-^ABS(V(L)  )+ABS(W(L)  )+AlX(L) ) 5803*5999,5803 

5803  1X2=1 
GO  TO  5999 

5999  IF ( K-K2) 4999#  5804  #5804 

5804  IF (KZ2) 4999# 5805# 4999 

5805  IF(ABS(U(L) ) +ABS ( V (L) ) +ABS ( W <L) ) +AIX ( L> ) 5806 # 4999, 5806 

5806  KZ2=1 
GO  TO  4999 

4999  IF(KZ1)5300#5222»5300 
5222  IF(K1-1)5300, 5300, 5000 

5000  IF(K1-K)5300, 5001, 5001 

5001  IF(ABS(U(L) )+ABS(V(L) )+ABS(W(L) )+A IX (L) ) 5002# 5300, 5002 

5002  KZ1=1 

5300  IF(JY1)5600, 5304*5600 
5304  IF(J1-1)5600, 5600, 5301 

5301  IF(J1-J)L600, 5302, 5302 

5302  IF(ABS<U(L) )+ABS(V(L) ) +ABS < W (L> ) +AI X (L> > 5303 » 5600 , 5303 

5303  JY1=1 

5600  IF(IX1)3342#5604#3342 
5604  IF ( I 1-1 ) 3342 , 3342, 5601 

5601  IF ( 1 1-1 ) 3342  #5602,5602 

5602  IF(AbS(D(L) )+A8S<V(L) ) +ABS <W (U ) +A I X (L> ) 5603 » 3342, 5603 

5603  1X1=1 
GO  TO  3342 

CAME  HERE  BECAUSE  THE  CELL  IN  QUESTION  (L)  IS  VOID 

*  3340  PRR=0. 

URR=U (L+l ) 

PABOVE-0  * 

VA0OVE-V (N) 

PZR=0. 


OOO  DO  DO  ODD  OO  OO  O  O  O  O  O  DO  OO 


nZR=W (NN) 

C 

C  SET  THL  AbOVE  QUANTITIES  TO  BELOW 

3342  VbL0=VAt30VE 
PULO=PABOVE 

SET  ThL  RIGHT  QUANTITIES  TO  THE  LEFT 

PL( J)=PRR 

UL(J)=URR 

SET  the  front  quantities  to  behind 

PUIND(M)=PZR 

UBIND(M)=WZR 

update  the  indices 

L=N 

M=  M+IMAX 

TERMINATION  OF  LOOP  ON  J - (Y) 

5  CONTINUE 

CHECK  ON  ADVANCING  OR  DECREASING  GRID  COUNTERS 

5700  LJ=L-IMAX 

IF ( JY2 ) 4  *  5701  *  4 

5701  IF(ABS<U(LJ) )+ABS(V(LJ) )+ABS(W(LJ) )+AIX(LJ) ) 5702* 4* 5702 

5702  JY2=1 
GO  TO  4 

TERMINATION  OF  LOOP  ON  I - (X ) 

4  CONTINUE 

termination  OF  LOOP  ON  K - <Z) 

3  CONTINUE 


CHECK  FOR  FIRSi  OR  SECOND  PASS 

44  IF(VEL-1.)46,45»46 

45  VEL=Q. 

RECYCLE 
GO  TO  2 

HAVE  COMPLETED  BOTH  PASSES 

46  CONTINUE 

increase  or  decrease  counters  as  requred 

11=11-1X1 

5900  1F(I1)5901*5901*5902 

5901  11=1 

5902  12=12+1X2 

5903  IF ( 1 2- 1 MAX) 5905 »  5905 » 5904 

5904  I 2= I MAX 

5905  Jl= Jl-UY 1 
IF(J1)59G6*5906*59Q7 

5906  Jl=l 

5907  J2=J2+JY2 

IF (J2-UMAX) 5909»  5909*  5908 


oooooooooooooo 
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5900  J2=JMAX 

5909  K1=K1-KZ1 
IF(K1)5910#5910»5911 

5910  Kl=l 

5911  K2=K2+KZ2 

IF (K2-KMAX) 5913# 5913» 5912 

5912  K2=KMAX 

5913  RETURN 
8999  NK=11 

GO  To  9999 

9902  NK=20 

GO  To  9999 

9900  NK=3305 

GO  To  9999 

9901  NK=330fa 

GO  To  9999 

9903  NK=3307 

GO  To  9999 
9911  NK=3316 

GO  TO  9999 

9904  NK=3310 

GO  TO  9999 

9905  NK=3320 

GO  TO  9999 

9906  NK=32 

GO  TO  9999 

9907  NK=332S 

GO  TO  9999 
9900  NK=4420 

GO  To  9999 

9909  NK=432*+ 

GO  TO  9999 

9910  NK=4326 
9999  NR=3 

WRITE (6 #8500) I# J#K#L»M»N#NN#NK»NR 
8500  FORMAT (916) 

CALL  UNCLE 
C  CALL  DUMP 
CALL  DUMP 
END 

(31  FOR  PH2/S»PH2/S»PH2/SS 

SUBROUTINE  PH2  PH2  QOIq 

DIMENSION  AIX (6000) »AMX (6000) »U (6000) #V (6000) #W(6000) »P(6000) » 

*********  3DOIL  *************************** 


HERE# WE  APPROXIMATE  THE  TRANSPORT  TERMS  LEFT  OUT  OF  THE 
MOMENTUM  AND  ENERGY  EQUATIONS  IN  PHI  BY  MOVING 
MASS# (SOLVING  THE  MASS  CONSERVATION  EQUATION)  THIS  MASS 
THEN  CARRIES  ENERGY  AND  MOMENTUM  ACROSS  THE  FIXED 
GRID  LINES 

AMPY  =  MASS  FLOW  AT  THE  TOP 

AMUT  =  X  MOMENTUM  COMPONENT  OF  THIS  MASS 

AMVT  =  Y  MOMENTUM  COMPONENT  OF  THIS  MASS 

AMWT  =  Z  MOMENTUM  COMPONENT  OF  THIS  MASS 

OELET=  SPECIFIC  EnERGY  OF  THIS  MASS 


oooo  on  ooo  non  oooooooooooooooonooooooononnnoo 


mMMP  =  MASS  FLOW  AT  THE  RIGHT 
AMUR  =  X  MOMtHTUM  COMPONENT  OF  THIS  MASS 

AMVR  =  Y  MOMENTUM  COMPONENT  OF  THIS  MASS 

AMWR  =  2  MOMENTUM  COMPONENT  OF  THIS  MASS 

UELER=  SPECIFIC  ENERGY  OF  THIS  MASS 

ammy=  mass  flow  at  the  bottom 

AMMUr  X  MOMENTUM  COMPONENT  OF  THIS  MASS 

AMMV=  Y  MOMENTUM  COMPONENT  OF  THIS  MASS 

AMMW=  2  MOMENTUM  COMPONENT  OF  THIS  MASS 

OELEb=  SPECIFIC  ENERGY  OF  THIS  MASS 

GAMC  =  MASS  FLOW  AT  THE  LEFT 
FLEFT=  X  MOMENTUM  COMPONENT  OF  THIS  MASS 
YAMC  =  Y  MOMENTUM  COMPONENT  OF  THIS  MASS 

2M0M  =  W  MOMENTUM  COMPONENT  OF  THIS  MASS 

sigc  =  specific  energy  of  this  mass 

BMASS=  MASS  FLOW  AT  THE  BACK 

BXMOM-  x  momentum  component  of  this  mass 

oYMOMs  Y  MOMENTUM  COMPONENT  OF  THIS  MASS 

d2M0M=  Z  MOMENTUM  COMPONENT  OF  THIS  MASS 

BENR  =  SPECIFIC  ENERGY  OF  THIS  MASS 

FMASS=  MASS  FLOW  IN  FRONT 
FXMOM=  X  MOMENTUM  COMPONENT  OF  THIS  MASS 

FYMOM=  Y  MOMENTUM  COMPONENT  OF  THIS  MASS 

FZM0M=  2  MOMENTUM  COMPONENT  OF  THIS  MASS 

FENR  =  SPECIFIC  ENERGY  OF  THIS  MASS 

RE2=0. 

INITIALI2E  THE  FLAGS  FOR  ADVANCING  THE 
ACTIVE  GRID  COUNTERS  TO  ZERO 
1X1  =  0 
1X2=0 
JY1  =  0 
JY2=0 
K2l  =  0 
K22=0 
SUM=0 . 

CALL  SLITE(O) 

UO  LOOP  ON  K 

5  DO  1  K=K1#K2 
LL=(K-1)*IXMAX 

DO  LOOP  ON  I 

6  uO  2  1=11 r 12 

NOTE  (L)  IS  THE  CELL  INDEX  =<J-1)IMAX  +  I 
+  (K-l)IXMAX  (NOTE  IXMAX=( IMAX) ( JMAX) 


C 


L=LL+I+(J1-1)*IMAX 


c 

c 


7 


DO  LOOP  ON  J 


{ 


r 


k 


c 

c 

C* 


DO  3  J=J1,J2 

NN  =  THL  INDEX  OF  THE  CELL  IN  FRONT  =  L+  IXMAX 

nn=l+ixmax 

n=  the  index  of  the  cell  above  =  l+  imax 
n=l+imax 


c 

c  m=  index  of  the  cell  in  question  for  a  single  plane 
c 

M=(J-1)»IMAX+I 

C  N1»N2»N3»N4»N5*NF  ARE  FLAGS  TO  SET  BOUNDARY  CONDITIONS 
C  AT  THE  6  FACES  OF  THIS  GRID 
C 

C  N1  REFERS  TO  THE  LEFT 

C  N2  REFERS  TO  THE  RIGHT 

C  N3  REFtRS  TO  THE  TOP 

C  N4  REFERS  TO  THE  BOTTOM 

C  N5  REFERS  TO  BEHInD 

C  No  REFERS  TO  IN  FRONT 

C 

C  FREE  SURFACES  WITHIN  THE  GRID  ARE  TREATED  AS  FOLLOWS  » 

C  -  if  the  mass  flow  into  a  empty  cell  produces  a 

C  DENSITY  THAT  IS  LESS  THAN  T020NE(A  INPUT  NUMBER  LIKE 

C  .001  OF  RHONOT) *  THE  MASS  FLUX  IS  SET  TO  ZERO 

C  . 

*  600  IF ( J-l) 9903 » 601 #9302 
9302  IF(J1-J)9#603#9 


C 

C  BOTTOM  BC.  HAS  BEEN  SET 
C 

601  IF ( AMX (L) ) 99 04*603*602 
C 

c  we  are  at  the  bottom  of  the  x-y  plane 

602  IFIV(L) )604»603»603 
C 

c  set  y  component  Of  momentum  to  0. 

603  AMMV=0. 

GO  TO  698 

C 

C  CALCULATE  THE  MASS  FLUX  AT  THE  BOTTOM 

604  AMMY=AMX(L)/DY(J)*V(L)*DT 
C 

C  CHECK  SO  MASS  FLUX  DOES  NOT  MORE  THAN  EMPTY  THE  CELL 

605  IF(AMMY+AMX(L) ) 9300  *  607, 607 
9300  AMMYr-AMX(L) 

607  IF(N4)6Q9, 608*609 


C- 

C*  BOTTOM  80UNDARY  IS  TRANS 

608  AMMU=AKMY*U(L) 

C: 

C  CALCULATE  THE  3  MOMENTASrTHE  ENERGY 

C  SUBTRACT  THIS  ENERGY  LOSS  FROM  ETH 
AMMV=AMMY*V(L) 

ammw=ammy»w(l) 

WS=(U(L)**2+V(L)**2+W(L)**2>/2. 


oooo  oooo  o  o  o  o  ooo  oooo 


□ELEB=AIX(L)+WS 

ETH=ETH+AMMY*DELEb 

IF(-AMMY/(DX(I)*DY( J)*DZ(K> )-Z(80) ) 610*610*6600 
'  6o00  REZ=1.0 

uO  TO  blO 

BOTTOM  BOUNDARY  Is  REFLECTIVE*  NET  MOMENTA  CHANGES  2MV 
o09  AMMV=2.*AMMY*V(L) 

SET  MASS *  X  AND  Z  MOMENTA  AND  SFECIFIC  ENERGY  TO  0. 

69d  AMMY=0. 

AMMU=0. 

AMMW=0« 
deleb=o. 

610  CONTINUE 

***  FINISHED  WITH  THE  BOTTOM  BC.  *********** 
****************************************************************** 
9  IF ( I-l ) 8999 *10*9301 
9301  IF(I1-I)5G6>9310»506 
9310  JJ=J 

GO  TO  20 

WE  ARE  ALONG  THE  LEFT  BOUNDARY  (  1=1) 

10  NL=L 

11  JJ=J 

15  IF ( AMX (NL) ) 9900  * 20*16 

20  FLEFT(UJ)=0. 

GO  TO  5504 

16  IF(U(NL) )17o20»20 

CALCULATE  MASS  FLUX 

17  GAMC ( JJ) =AMX (NL) /DX ( 1 ) *U (NL) *DT 

21  IF(GAMC(JJ)+AMX(NL) )22*500*500 

22  GAMC ( JJ)s-AMX (NL) 

500  IF(N1)501»502*501 

LEFT  BOUNDARY  IS  TUANS 
CALCULATE  THE  3  MOMENTAS*THE  ENERGY 
SUBTRACT  THIS  ENERGY  LOSS  FROM  eth 
502  FLEFT(JJ)=GAMC(JJ)*U(NL) 

YAMC ( JJ) =GAMC ( JJ) *V(NL) 

ZMOM(JJ)=GAMC(JJ)*W(NL) 

W s= ( U ( NL ) **2+ V ( NL ) **2+W (NL)**2)/2. 

SIGC ( JJ) =AIX (NL) +WS 
ETH=ETH+GAMC ( JJ) *SIGC ( JJ) 

IF (“GAMC (J)/(DX(I)*DY(J)*DZ(K) ) -Z (75) ) 503* 503*6610 
6610  REZ=1.0 
GO  TO  503 

LEFT  BOUNDARY  IS  REFLECTIVE*  NET  MOMENTA  CHANGE  =2MU 

501  FLEFT(JJ)=2.*GAMC(JJ) *U(NL) 

SET  MASS*  Y  AND  Z  MOMENTA  AND  SPECIFIC  ENERGY  TO  0* 

5504  GAMC(JJ)=0. 

YAMC ( JJ ) =0 • 

ZMOM ( JJ) =0 • 


o  o  o  o  o  .  p  n  oo  nnn  on  oo  o-r>  o  <o  on  on  o  o  o  o  o 


siec(ju)=o. 

503  CONTINUE 
103  CONTINUE 

FINISHED  WITH  LEFT  BOUNDARY  CONDITIONS 

***  *  *  41  *  **  *  *  *  *  *  *  **  *  *  *  *  *  *  *  *  *  4c  *  *  41 41  *  *  *  *  *  *  *  *  *  *  *  #  *  *  #  *  #  *  *  *  *  *  *  #  *  *  *  *  *  *  *  ♦  *  *  * 

506  IF<K-l)990l#23»9303 
9303  IF(K1-K)31#250#31 

23  IF(Amx(L) )9902»25o»24 
250  BZMOM(M)=0. 

GO  TO  25 

CHECK  2  COMPONENT  OF  VELOCITY 

24  IF ( W (l) ) 26 » 250 » 250 

SET  THE  5  DATA  BEHIND  TO  0. 

25  BMASS(M)=0. 
bXMOM(M)=0. 

UYMOM(M)=0. 

BENR(M)=0. 

GO  TO  31 

VELOCITY  IS  -  f CALCULATE  THE  MASS  FLUX 

26  DMASS(M)=AMX(L)/02(K)*W(L)*DT 

CHECK  SO  WE  DONT  EMPTY  MORE  MASS  THAN  THERE  IS 

27  IF(BMASS(M)+AMX(L) )28#40»40 

28  BMASS(M)=-AMX(L) 

CHECK  FOR  TRANS  OR  REFLECT 

40  IF (N5)41#29»41 

reflective 

41  BZM0M(M)=2.*BMASS(M)*W<L) 

GO  TO  25 

TRANSMITTIVE 

CALCULATE  THE  MOMENTAS  OF  THIS  MASS 

29  BXMOM(M)=BMASS(M)*U(L) 

BYMOM(M)=BMASS(M)*V(L) 

UZMOM(M)=BMASS(M)*W(L) 

CALCULATE  THE  TOTAL  ENERGY  CARRIED  BY  THIS  MASS 

30  WS=<U(L)**2+V(L)**2+W(L)**2)/2. 

UENR(M)=AIX(L)+WS 

REMOVE  THE  ENERGY  LOSS  FROM  ETH 
ETH=ETH+BMASS(M) *BENR (M) 

IF(-BMASS(M)/(DX(I)*DY(J)*DZ(K) )-Z(78) ) 31 r 31,6620 
6620  REZ=1.0 

HAVE  CALCULATED  THE  DATA  BEHIND#  NOW  CHECK  ON  JMAX 

***>**♦*******************************************♦** 

NOW#  UP  TO  THIS  POINT#  WE  HAVE  TAKEN  CARE  OF 


ooooooo  oo  non  •  r>  o  no  oo  oo  non  noon  r/ non 
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BOTH  REFLECTIVE  AND  TRANSMITTIVE  BOUNDARIES 
AT  Th£  BOTTOM*  LEFT *  AND  BACK 
***************************************************** 

* ********** **************************************** *************** 


VABOVE  CALC. 

********************************** ********** 
31  VEL=U. 


set  up  to  calculate  vabove 

IF ( JMAX-J) 211  *  211 *  207 

WE  are  at  the  top  of  the  grid 

211  VEL=1. 

GO  To  208 

check  ^ell  above 

207  IF(AMXtN) )215*215,2l4 

CELL  ABOVE  IS  VOID 

214  IFtAMX(L) >216*216,209 

CELL  (L)  IS  VOID*  BUT  CELL  ABOVE  IS  OCCUPIED 
216  VAB0VE=V(N) 

GO  TO  212 

215  IF(AMXU-)  >205*205,208 

BOTH  CELLS  ARE  VOID 

205  VABOVE=0. 

GO  To  212 

208  VA80VE=V(L) 

GO  TO  212 

BOTH  CELLS  ARE  OCCUPIED 

209  VABOVE=(V(L)+V(N) )/2. 

212  FS=0. 


U  RIGHT  CALC. 

******************************************** 

NOW* BEGIN  calculation  of  urr 

404  IF (IMAX-I) 412*412, 405 

405  IF(AMX(L+1) )411*411#409 

409  IF(AMX<L> >410,410*407 

C  CELL  (L)  IS  VOID  *BUT  CELL  TO  THE  RIGHT  IS  FILLED 

410  URR=U(L+1) 

GO  TO  408 

C 

C  MASS  ON  THE  RIGHT=0. 

411  IF(AMX(L> >403*403,406 
403  URR=0. 
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60  TO  40S 
C 

C  WE  ARE  AT  THE  RIGHT  OF  THE  GRID 

'  412  FS=1. 

406  URR=U(L) 

GO  TO  403 

.  C 

*C  THIS  IS  THE  NORMAL  PATH 

C  BOTH  CELLS  ARE  FILLED 

407  URR=(U(L)+U(L+l))/2. 

408  CONTINUE 
C 

C  W  IN  FRONT  CALC. 

C  ******************************************** 

c 

C  NOW  ,WE  HAVE  VABOVE  (V  AT  THE  TOP)  AND  URRtV  AT  THE  RIGHT) 

5503  AREA=0. 

C 

C  LETS  CALCULATE  WOUT  (THE  Z  COMPONENT) 

504  IF(KMAX-K)512»512»505 

505  IF(AMX(NN) ) 51 lr 51 lr 509 

509  IF(AMX(L) )5l0,510,507 
C 

.C  CELL  (L)  IS  VOID  ,BUT  CELL  IN  FRONT  IS  OCCUPIED 

510  WOUT=W(NN) 

GO  TO  503 

.C 

C  MASS  IN  FRONT  IS  0.  .  \ 

511  IF(AMX(L))5513, 5513, 5066 
5513  WOUT=0. 

GO  TO  508 
C 

C  WE  ARE  AT  THE  FRONT  BOUNDARY 

512  AREA=1. 

5066  WOUT=W(L) 

GO  TO  508 
C 

C  NORMAL  PATH 

C  BOTH  CELLS  ARE  FILLED 

507  WOUT=(W(L)+W(NN) )/2. 

C 

C  IF  VEL  IS  GREATER  THAN  0.  »WE  ARE  IN  THE  TOP  CELL  <J=GMAX) 

C  IF  FS  IS  GREATER  THAN  0.,  WE  ARE  IN  THE  RIGHT  COLUMN 

C  THAT  IS  (I=IMAX) 

C  IF  AREA  IS  GREATER  THAN  0.,  WE  ARE  IN  THE  FRONT  SLAB(X-YPLANE) 

C  THAT  IS  (K=KMAX) 

508  CONTINUE 
C 

c  *************************************************** 

C  NOW*  FINALLY,  we  HAVE  ALL  3  INTERFACE  VELOCITIES 
C  ************** 

C  HERE  WE  BEGIN  THE  CALCULATION  OF  THE  3  FLUXES 

C 

C  ******************************************** 

C  IF  TOP  IS  TRANS.  THE  ENERGY  LOSS  IS  CALCULATED  LATER 

C  ******************************************** 

100  IF(VAB0VE)102, 101, 1103 


c 

C  Y  FLUX  IS  POSITIVE 

1103  IF(AMX(L) ) 9910 » 101 ,104 
•  C 

c  SET  INDICES 

104  IF (J-l ) 9910 » 6104 » 6105 
-‘6105  KP=L-IMAX 

IF(Amx(KP) )9910»6l06»6104 
6106  IF ( ABS ( VABOVE-VELOC ) /VELOC-BUG ) 1 07 » 61 04 , 6104 
6104  LY=L 
JY=J 

IF (VEL) 105 » 105 » 109 
C 

C  Y  FLUX  =  0. 

101  AMPYrO. 

AMUT=0. 

AMVT=0. 

AMWT=0« 

DELET=0. 

GO  TO  115 

C 

C  Y  FLUX  IS  NEGATIVE  (DOWN) 

C  FLUX  IS  NEGATIVE  FROM  CELL  (L) 

102  IF (VEL) 106*106*101 
-C 

C  FLUX  IS  -  BUT  TROM  THE  TOP  CELL 
C  CHECK  CELL  ABOVE 

-  106  IF(AMX<N) )9911,10l,107  \ 

C 

C  FLUX  IS  NEGATIVE»BUT  CELL  MASS  =0, 

107  LY=N 
JY=J+1 

105  VABOVE=  (VABOVE)/(l.+(V(N)-V(L) )/DY(JY)*DT/SBOUND) 

C 

C  CALCULATE  FLUX  AT  THE  TOP 

C  *********** 

109  AMPY=AMX (LY) *VABOVE/OY ( JY) *OT 

110  IF (VEL) 115 » 115 » 111 
C 

C  WE  ARE  AT  THE  TOP, CHECK  THE  BC.  AT  THE  TOP 

111  IF (N3) 112, 115* 112 
C 

C  REFLECT  THE  MASS 

112  IF ( AMPY) 115» 115» 113 

113  AMVT=-2,*AMPY*V(L) 

114  AMPY=0. 

AMUT=0. 

AMWTzO • 

DELET=0. 

C  ******************************************** 

C  IF  RIGHT  IS  TRANS.  THE  ENERGY  LOSS  IS  CALCULATED  LATER. 

.C  ******************************************** 

C 

C  BEGIN  CALCULATING  THE  FLUX  AT  THE  RIGHT 

115  IF (URR) 118'116'117 

116  AMMP=0. 

AMUR=0 • 


oooooo  o  o  o  o  oo  non  oooo  o  r>  ’  o  o  o 


AMVR=0. 

AMWR-0. 

DELERsO. 


CALCULATE  THE  MASSFLOW  1^  THE  Z  DIRECTION 
60  TO  300 

117  IF(AMX<L) }99llrll6’120 

120  LX=L 
IX=I 

IF (FS) 119» 119» 122 

RIGHT  FLUX  IS  NEGATIVE 

118  IF(FS) 130 » 130 » 116 

FLUX  IS  NEGATIVE  »BUT  CELL  MASS  IS  ZERO  ALSO 

CHECK  THE  CELL  TO  THE  right 
130  IF(AMX(L+1) )9911#116#121 

121  LX=L+1 
IX=I+1 

119  WS=(U(L+1)-U(L) ) 7DX ( IX) *DT/SB0UND 
URR=URR/(1.+WS) 

CALCULATE  THE  MASS  FLUX  AT  THE  RIGHT 
***************** 

122  AMMP=AMX(LX)/DX(IX)*URR*DT 

IF(FS) 123»300* 123  \ 

CHECK  the  BOUNDARY  CONDITION 

123  IF<N2)124,300»124 

Trans 

124  IF  C AMMP) 300r  300* 125 

reflect  the  MASS 

125  AMUR=-2.*AMMP*U(L) 

AMMP=0* 

AMVR=0. 

AMWR=0. 

DELER=0. 

IF  FRONT  IS  TRANS.  THE  ENERGY  LOSS  IS  TAKEN  CARE  OF  LATER. 

DO  THE  Z  COMPONENT  NOW 

SET  THE  5  OATA  IN  FRONT  TO  0. 

300  IF(W0UT)318»316»317 

316  FMASS=0. 

FXMOM=0. 

FYMOM=0. 

FZMOM=0. 

FENR=Q» 

GO  TO  700 

317  IF(AMX(L) )9912»316»3200 
3200  LZ=L 

IZZ=K 


o  o  o  n  o  o  on  oooo  non  o  o  ooo 


IF (AK£A) 319 #319 #322 
C 

C  FRONT  FLUX  IS  NEGATIVE 

318  IF ( AREA) 320# 320 #316 

FLUX  IS  NEGATIVE#  BUT  FROM  IN  FRONT 
CHECK  CELL  IN  front 

320  IF (AMX(NN) ) 9912# 316# 321 

FLUX  IS  NEGATIVE  #BUT  CELL  MASS  =0. 

321  LZ=NN 
IZZ=K+1 

319  WS=(W(NN)-W(L) )/DZ( IZZ) *DT/SB0UND 
W0UT=W0UT/(1.+WS> 

CALCULATE  THE  MASS  FLUX  IN  FRONT 
*****♦***♦*>11*  *  * 

322  FMASS=AMX (LZ) /DZ ( IZZ) *WOUT*DT 
IF (AREA) 323 #700  #  323 

323  IF(N6) 324 #700#  324 

TRANS 

REFLECT# WE  ARE  IN  FRONT 

324  IF {FMASS) 700»700  #  325 

REFLECT  THE  MASS 

325  FZM0M=-2.*FMASS*W(L) 

FMASS=0. 

FXMOM=0. 

FYMOM=0. 

FENR=0» 


NOW  WE  HAVE  ALL  3  FLUXES  AND  ALL 
THE  BOUNDARY  CONDITIONS  HAVE  BEEN  SET 
700  IF(A,v(pY)760#980#76l 

TOP  FLUX  IS  - 

760  IF(AMPY+AMX(N))762#980#980 

762  AMPYr-AMX(N) 

GO  TO  980 

761  IF (-AMPY+AMX (L) ) 763#  980#  980 

763  AMPY=AMX(L> 

980  IF(AMMP)7300»981#7301 

7300  IF(AMmP+AMX(L+1) )7302#981#981 

7302  AMMP=-AMX(L+1) 

GO  TO  981 

7301  IF(-AMMP+AMX(L) ) 7303# 981 # 981 

7303  AMMP=AMX(L> 

981  IF (FMASS >7400 #982, 7401 

7400  IF(FMASS+AMX(NN) )7402#982#982 

7402  FMASS="AMX(NN) 

GO  TO  982 

7401  IF (“FMASS+AMX (L) ) 7403 #982 #982 

7403  FMASS=AMX(L) 

982  WS=GAMC(J) 


W '  1  '*er 


IF( WSJ  902 #901 #901 
901  WS=0. 

.902  WSA=AMMY 

IF(WSA)904#903»903 
903  WSA=0. 

,904  WSB=BMASS(M) 

IF(WSB)906#905#905 
905  WSB=0. 


906  WSC=AMX(L)+WS+WSA+WSB 

907  WS=AMPY 

IF (WS) 950 #950  #909 

950  WS1=0. 

IF(K-K1)951»953#951 

951  IF(GamC<J+1))952#953»953 

952  WS1=GAMC(J+1) 

953  WS2=0. 

NA=M+IMAX 

IF ( K-Kl ) 955  #  957  r  955 

955  IFCBMASS(NA) )956»957#957 

956  WS2=GMASS(NA) 

957  WS3=WS1+WS2+AMX(N) 

958  IF(AMPY+WS3)959,908»908 

959  AMPY=-WS3 
GO  TO  908 

*  908  WS=0. 

909  WSA=AMMP 

,  .  IF(WSA)970»970»9ll 

*  970  WS1=0. 

IF(K-K1)971#973#971 

971  IF<BMASS(M+1))972,973#973 

972  WS1=BMASS(M+1) 

973  WS3=WSl+AMX(L+l) 

974  IF(AMMP+WS3)975»910#910 

975  AMMP=-WS3 
GO  TO  910 

910  WSA=0. 

911  WSB=FMASS 
IF(WSB)912»912»913 

912  WSB=0. 

913  WST=WS+WSA+WSB 
IF(WST)921»921»931 

931  IF(WSC)932#932»933 

932  IF ( AMPY ) 934 » 934 » 935 

935  AMPY=0. 

934  IF ( AMMP ) 936 # 936 » 937 

937  AMMPzQ • 

936  IF(FMASS)921»921»938. 

938  FMASS=0t 

.  GO  TO  921 

*  933  IF(WSC-WST)914#921#921 

914  WSU=WSC/WST 

.  WS=WS*WSD 

WSA=WSA*WSD 

WSB=WSB*WSD 

915  IF (WS) 917 #917 #916 

916  AMPY=WS 

917  IF (WSA) 919 #919 #918 


onn  n  o  on  on  no 


918  AMMP=WSA 

919  IF ( WSB) 921 » 921 #  920 

920  FMASS=WSB 
'.921  CONTINUE 

IF  (  ANiPY  ) 703#  2700 »  /  02 
2700  IF ( J-JMAX) 717 #  2701 » 2701 
§701  IF (N3) 716 #717 #716 
717  AMUTsO# 

AMVTzO • 

AMWT=0. 

uelet=o. 

GO  TO  716 
C 

C  TOP  FLUX  IS  + 

702  IF(JF,AX-J)9913#704#705 

CHECK  CELL  ABOVE 

705  IF(Af-.X(N)  )9914#706»704 

FREE  SURFACE  AT  TOP 

706  IF ( AMP Y/ (UX ( I ) *D Y ( J ) *DZ (K ) ) -TOZONE >  707 » 704 » 704 

DENSITY  IS  TOO  SMALL#  SET  FLUX  =0. 

•707  AMPYsO. 

GO  To  717 

TOP  FLUX  IS  NEGATIVE 

703  IF(J-1)9913»701»709 

709  IF(AMX(L) )9914#710»701 

710  IF(-AMPY/(DX(I)*DY(J)*DZ(K))-TOZOnE)711#701#701 

711  AMPYsO. 

GO  TO  717 

ADD  UP  THE  MASSES  (REMEMBER  THEY  HAVE  DIRECTION) 

704.  DELMzgAMC ( J>  +AMMY+BMASS(M) -AMPY 
IF(VEL)9914#712#720 
C  AT  TOP  OF  GRID 
720  IF(N3)713#714#713 
C 

C  TOP  boundary  is  transmittive 

714  JS-  U (L ) **2+V (L) **2+W (L) **2 
ETH=ETH-AMP Y* ( A I X ( L ) + WS/2 , ) 

IF(AMPY/(DX(I)*DY(J)*DZ(K) )-Z(79) )712#712#6630 
6630  R£Z=1. 

C 

C  CALCULATE  THE  MOMENTUMS 

712  AMUT=AMPY*U(L) 

AMVT=AMPY*V(L) 

AMWT=AMPY*W(L) 

GO  To  713 

:  701  AMUT=AMPY*U(N) 

AMVT=AMPY*V(N) 

AHWT=AMPY*W(N) 

0cLET=AIX(N)+(U(N)**2+V(N)**2+W(N)**2)/2. 

716  D£LM=GAMC ( J ) +AMM Y+BMASS ( M ) -AMPY 
GO  TO  715 


on  no  o  o  on  on  o  on  on  o  oooorv 


© 


713  0ELET=AIX(L)+(U(L)**2+V(L)**2+W<L>**2)/2. 

715  SIGMU=FLEFT ( J) +AMMU-AMUT+BXMOM ( M> 
SIGMV=YAMC(J)+AMMV-AMVT+BYMOM(M) 

DELEK=G AMC ( J ) *S I GC ( J ) +AMMY*DELEB-AMP Y*DELET*BMASS ( M ) *BENR ( M ) 
S I GMW=ZMOM ( J ) +AMMW-AMWT+BZMOM ( M  > 

GO  TO  7000 

NOW » DO  THE  SAME  FOR  THE  X  DIRECTION 
MASS  FLUX  AT  THE  RIGHT 


7000  IF(AMMP) 7003*2702#  7002 

2702  IF( I- IMAX) 7017* 2703 »2703 

2703  IF(N2)7016'7017»7016 

FLUX  IS  POSITIVE 

7002  IF (  IMAX-I >9914*7004 #7005 
7017  AMUR=0. 

AMVR=0. 

AMWR=0. 

DELER=0. 

GO  TO  7100 

CHECK  CELL  TO  THE  RIGHT 

7005  IF(AMX(L+1) )9915*7006#7004 
FREE  SURFACE  AT  THE  RIGHT 

7006  IF ( AMMP/ (UX( I)*DY ( J)*DZ(K) > -T0Z0NE ) 7007 » 7004 » 7004 

OENSITY  IS  TOO  SMALL#  SET  FLUX  TO  0. 

7007  AMMP=0 • 

GO  TO  7017 

RIGHT  FLUX  IS  NEGATIVE 

7003  IF(I-1)9914#7001»7009 

7009  IFCAMX(L) )9915»7010#7001 

7010  IF ( -AMMP/ ( DX ( I ) *DY ( J ) *DZ ( K ) ) -T0Z0NE ) 701 1 #  7001 # 7001 

FLUX  IS  TOO  SMALL 

7011  AMMP=0. 

GO  TO  7017 

7004  OELMsDELM-AMMP 
IF (FS) 9915 #7012 #7020 

AT  RIGHT  BOUNDARY 
7020  IF(N2)7l50»7014»7l50 

RIGHT  BOUNDARY  IS  TRANSMITTIVE 
.  7014  WS=U(l)**2+V(L)**2+W(L)**2 
ETH=ETH-AMMP* ( A I X ( L ) +WS/2 . ) 

IF(AMMP/(DX(I)*DYCU)*DZ(K)-Z(77))7012#7012»6640 
.  6640  REZ=1, 

*  7012  AMUR=AMMP*U<L) 

AMVR=AMMP*V(L) 

AMWR=AMMP*W(L) 

DELER=AIX(L)+(U(L)**2+V(L)**2+W(L>**2)/2. 

GO  TO  7150 


7001  AMUfv=AMMP*U(L+l) 

AMVR=aMiMP*V(L+1> 

AMWR=aMMP*W(L+1> 

DELER=A IX(L+1)  +  (U(L+1) **2+ V <L+ 1  >  **2+W ( L+l ) **2) /2. 
7016  DELM=DELM-AMMP 
C 

•C  SUM  UP  TOTAL  MOMENTA 
7150  SIGMU=SIGMU-AMUR 
SIGMV=SIGMV-AMVR 
SIGMW=SIGMW-AMWR 

delek=delek-ammp*deler 

C 

c  do  the  same  for  the  flux  in  front 

C  "HI.  Mini  Ml.  If  . . .  till  Mttlll 

c 

7100  IF(FMASS) 7103 #2704 #9982 

2704  IF (K-KMAX) 7117 #2705 #2705 

2705  IF(N6)4000#7117#4000 
C 

C  FLUX  IS  POSITIVE 

9982  IF(KMAX-K)9916#7104»7105 
7105  IFCAMXCM.N))9917»7106#7104 
C 

C  FREE  SURFACE  IN  FRONT 

*  7106  IF(FMASS/(DX(I)*DY(J)*D2(K) )-TOZONE)7l07»7104#7104 
C 

.C  DENSITY  IS  TOO  SMALL#  SET  MASS  =0. 

*  7107  FMASS=0.  \ 

7117  FXMOM=0. 

FYMOM=0 . 

fzmom=o« 

FENR=0. 

GO  TO  4000 
C 

C  FRONT  FLUX  IS  - 

7103  IF(K-1>9916#7101»7109 

7109  IF(AMXCL) )9917»7110»7101 

7110  IF(-FMASS/(DX(I)*DY(J)*DZ(K) )-T0Z0NE) 7111 #7101 #7101 
C 

C  FLUX  IS  TOO  SMALL 

7111  FMASS=0. 

GO  TO  7117 

7104  DELM=DELM-FMASS+AMX(L) 

IF (AREA) 9916»7112»7120 

7120  IF(N6 >8000/71 14 #8000 
C 

C  FRONT  BOUNDARY  IS  TRANSmITTIVE 

7114  WS=U(L)**2+V(L)**2+W(L)**2 
ETH=ETH-FMASS* ( A IX ( L ) +WS/2 . > 

IF ( FMASS/ ( DX ( I ) *DY ( J>  *DZ ( K ) ) -Z ( 76 >  >  71 12 # 7112 # 6650 
6o50  R£Z=1. 

7112  FXMOM=FMASS*U(L> 

FYMOM=FMASS*V(L) 

FZMOM=FMASS*W(L> 

FENR=aIX (L)  +  (U(L> ♦*2+V(L>  **2+W(L)**2)/2» 

GO  TO  8000 

7101  FXM0M=FMASS*U(NN) 


FYMOM=FMASS*V(NN> 

FZM0m=FMASS*W(NN) 

71.16  F£NR=A  I X  <  NN )  +  <  U  <  NN  >  **2+V  <  NN ) * *2+W  <  NN  >  **2  > /2 . 

4.000  DELM=DELM-FMASS+AMX(L) 

8000  sigmu=sigmu-fxmom 
sigmv=sigmv-fymom 

•  ‘  SIGMW=SIGMW-FZMOM 

oelek=delek-fmass*fenr 

c 

C  TOTAL  MASS  AT  CYCLE  N+l 
C  ********** 

IF(DElM)544»545»540 

544  IF ( AmX (L) *1 .E-6+DELM) 9918 » 545 » 545 

545  DELM=0. 

GO  TO  550 

540  WS=U(L>**2+V<L)**2+W(L>**2 
WS=WS/2 • 

EnK=AMX (L)*(AIX(L)+WS) +DELEK 
GO  TO  541 
C 

C  HERE  ,WE  CALCULATE  THE  3  NEW  CELL  (L)  VELOCITIES  BY 
C  CONSERVING  MOMENTUM'  WE  ALSO  CONSERVE  THE  TOTAL 

C  ENERGY  t  THE  TOTAL  ENERGY  LESS  THE  KINETIC  IS 

C  .  THAN  THE  NEW  SPECIFIC  internal  ENERGY. 

C* 

C  NEW  X  VEL  COMPONENT 

541  U(L)=(SIGMU+AMX(L)*U(L))/DELM 

d-  \ 

C  NEW  Y  VEL  COMPONENT. 

546  V ( L ) r ( S I GMV+AMX (L) *V (L) ) /DELM 
C 

C  NEW  Z  VEL  COMPONENT 

547  W(L)=(SIGMW+AMX(L)*W(L) )/DELM 

548  WS=U(L)**2+V(L)**2+W(L)**2 
C 

C  NEW  SPECIFIC  INTERNAL  ENERGY 

549  AlX<L)=ENK/0ELM-WS/2. 

IF(AIX<L)-TMASS)7500#7500»7510 

7500  SUM=SUM+AIX<L)*OELM 
AIX(L)=0. 

7510  IF<ABS<U(L) )-XMAX)7501»750l»7502 

7501  WS=U(L)**2 
SUM=SUM+DELM*WS/2. 

U(L)=0. 

7502  IF(ABS<V(L))-XMAX)7503#7503#7504 

7503  WS=V(L)**2 
SUM=SUM+DELM*WS/2. 

V(L)=0. 

7504  IF ( ABS t W (L) J-XMAX) 7505 » 7505 » 7506 

7505  WS=W(L)**2 
SUM=SUM+DELM*WS/2. 

.  W(L)=0. 

7506  IFCAIX(L) )4001#550»550 
4001  SUM=SUM+AIX(L)*OELM 

AIX(L)=0, 


550  AMX(L)=OELM 
5800  IF(I-I2)5999'5801'580l 


o  r>  o  o  o 
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5801  IF 1 1X2) 5999,5802,  5999 

5802  IF ( AbS (U ( L) ) +ABS ( V (L) )+ABS(W (L)>+AIX(L))5803» 5999,5803 

5803  1X2=1 

GO  TO  5999 

5999  IF(K-K2)4999, 5804, 5804 
.5804  IF (K 22) 4999, 5805, 4999 

*5805  IF ( AbS (U (L) )+ABS(V(L) )+ABS(W(L))+AlX(L>  >5806,4999,5806 
5806  KZ2=1 

GO  TO  4999 

4999  IF(KZ1)5300, 5222, 5300 
5222  IF(K1-1)5300, 5300, 5000 

5000  IF(K1-K)5300, 5001, 5001 

5001  IF(ABS(U(L))+ABS(V(L.))+ABS(W(L>  >+AlX(L>> 5002, 5300, 5002 

5002  KZ1=1 

5300  IF(JY1)5600, 5304, 5600 
5304  IF( J1-1J5600, 5600, 5301 

5301  IF( J1-JJ5600, 5302, 5302 

5302  IF ( AbS (U (L) )+ABS ( V (L) )+ABS(W(L))+AlX(L) >5303,5600,5303 

5303  JY1=1 

5b00  IF ( 1X1)3342,5604,3342 
5604  IF( 11-1)3342,3342,5601 

5601  IF(I1-I)3342, 5602, 5602 

5602  IF ( AbS CU (L) )+ABS(V(L) )+ABS(W(L>>+AlX(L>  >5603,3342,5603 ' 

5603  1X1=1 
*3342  CONTINUE 

551  IF(AMX(L) >9919,553,9980 

9980  IF ( AMX ( L ) / (OX ( I ) *DY ( J ) *DZ ( K >  > -TOZONE > 9981 » 552,552 
•9981  AMLOST=AMLOST+AMX(L)  \ 

WS=U(L>**2+V(L)**2+W(L)**2 
WSR=AMX(L)*(AIX(L>+WS/2.> 

SUM=SUM+WSR 

elost=elost+wsr 

XMLOST=XMLOST+AMX (L>  *U (L> 

YMLOST=YMLOST+AMX(L>*V(L> 

zmlost=zmlost+amx ( L ) *W ( l  > 

AMX(L)=0. 

553  AIX(L)=0. 

U(L>=0. 

V (L) =0 • 

W(L>=0. 

P(L>=0. 

HERE  THE  FLUX  DATA  FROM  THE  RIGHT  IS  SET  TO  THE  LEFT 

552  GAMC(J)=AMMP 
FLEFT(J>=AMUR 
YAMC(J>=AMVR 
ZM0M(J)=AMWR 
SIGC ( J) =OELER 

HERE  THE  FLUX  DAT A  FROM  THE  TOP  IS  SET  TO  THE  BOTTOM 

554  AMMY=AMPY 
AMMU=AMUT 

ammv=amvt 

AMMW=AMWT 

oeleb=delet 

C 


c 


HERE  THE  FLUX  DATA  FROM  IN  FRONT  IS  SET  TO  THE  BACK 

555  BMASS(M)=FMASS 
UYMOM(M)=FYMOM 
OXMOM(M)=FXMOM 
UZMOM(M)=FZMOM 
U£NR(M>=  FENR 

continue 
l=l+imax 

C  4SS$SSS$$$$$$  END  OF  J  LOOP  $$$$$$$$$$$$$S$$S$S$S$$$$$$$SSS 

3  CONTINUE 

5700  LJ=L-IMAX 
IF(JY2)4#5701#4 

5701  IF ( AbS (U (LJ) ) +ABS ( V (LJ) ) +ABS( W (LJ) ) +AIX ( LJ) ) 5702#  4 #5702 

5702  JY2=1 

4  CONTINUE 

556  CONTINUE 

C  END  OF  I  LOOP  $$$$$$$$$$$$$$$$$$ 

2  CONTINUE 

557  CONTINUE 

c  end  of  k  loop  $$$$$ss$$$$$$$$$$ssss$$$ss$ 

i  continue 

8001  ETH=ETH-SUM 
ENE6=ENEG-SUM 
11=11-1X1 

*  5900  IF(I1)5901#5901#5902 
5901  11=1 

.5902  12=12+1X2 

•  5903  IF(I2-IMAX)5905#5905#5904  \ 

5904  I2=IMAX 

5905  J1=J1-JY1 

IF  (Jl) 5906 ’5906 ?590  7 

5906  Jl=l 

5907  J2=J2+JY2 

IF ( J2-UMAX ) 5909  #  5909  #  5908 

5908  J2=JMAX 

5909  K1=K1-KZ1 

IF  (Kl) 5910*5910  #  5911 

5910  Kl=l 

5911  K2=K2+KZ2 

IF (K2-KMAX) 5913# 5913# 5912 

5912  K2=KMAX 

5913  IF tREZ) 9950 #9950 #9951 
9951  CALL  REZONE 

9950  RETURN 

9903  NK=60Q 

GO  TO  9999 

9904  NK=601 

GO  TO  9999 
.  8999  NK=9 

GO  TO  9999 

9900  NK=15 

GO  TO  9999 

9901  NK=506 

GO  TO  9999 

9902  NK=23 

GO  TO  9999 
9910  NK=1103 


o  o  o 


GO  TO  9999 
9911  NK=l06 

GO  TO  9999 
'9912  NK=317 

GO  To  9999 

9913  NK=702 

GO  TO  9999 

9914  NK=705 

GO  TO  9999 

9915  i\IK=7005 

GO  TO  9999 

9916  NK=9982 

GO  TO  9999 

9917  NK=7l05 

GO  TO  9999 

9918  NK=544 

GO  TO  9999 

9919  NK=551 
9999  NR=4 

WRITE (6» 8500} I»J»K»L»M»N#NN»NK»NR 

WRITE (6#  8501 } GAMC ( J) #FLEFT ( J) » YAMC ( J) rZMOM(J) rSIGC<J) 

WR ITE ( 6 » 850 1 } AMMP » AMUR  #  AMVR » AMWR » DELER 
WRITE ( 6 , 8501 } AMMY , AMMU  #  AMMV » AMMW » DELEB 
WRITE ( 6  f 8501 } AMPY , AMUT » AM VT » AMWT » DELET 

WRITe(6»8501) BMASS (M) »BXMOM(M) #BYmOM(M) rBZMOM(M) »BENR<M> 
WRITE (6f 8501 )FMASS»FXM0M#FYM0M»F2M0M»FENR 
WRITE (6, 8501) AMX(L)»AIX(L) »U(L) »V (L) , W(L) rP(L) »CYCLE 
WRITE (6» 8501 )AMX(L+1)»U(L+1) r V<L+l) »W(L+1) 

WRITE (6» 8501 )AMX(NN)»U(NN) #V(NN) »W  CNN) 
WRITE(6»8501)AMX(L-1) »U(L-1) rV(L-l) rW(L-l) 

WRITE (6» 8501) AMX(N)»U(N) »V(N>  »W(N) 

LL=L-IMAX 

LBJ=L-IXMAX 

WRITE (6»  8501) AMX (LL) »U(LL) »V(LL) »W(LL) 

WRITE(6»8501) AMX(LBJ) »U(LBU) »V<LBJ> »W(LBJ> 

8501  FORMAT ( 1P8E12.5) 

8500  FORMAT (916) 

CALL  UNCLE 
CALL  DUMP 
END 

01  FOR  REZONE/S»REZONE/S»REZONE/SS 
SUBROUTINE  REZONE 

CHANGE  all  cell  dimensions  by  a  factor  OF  2. 

NOT£»  8  CELLS  BECOME  ONE  IN  THE  NEW  GRID 
CALCULATE  NEW  INOICES 
KKMAX=KMAX/2 
IIMAX=lMAX/2 
JJMAX=UMAX/2 
C  SET  UP  DO  LOOP  FOR  NEW  STORAGE 
KN1=-IXMAX 
DO  21  KKK=1»KMAX»2 
LL=(KKK-1)*IXMAX 
KN1=KN1+IXMAX 
KN=KN1 

DO  21  II=1»IMAX»2 

I=LL+II 

KN=KN+1 


o  o  o  o 


. 
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KK=KN 

DO  21  JJ=1 » JMAX 1 2 

k=i+ixmax 

J=I+IMAX 

l=k+imax 

WSA=AMX ( I ) +AMX ( 1+1 ) +AMX ( J ) +AMX ( J+l ) + 

1AMX(K)+AMX(K+1 ) +AMX (L)+AMX(L+1 ) 

IF<WSA>6*6»3 

C  CALCULATE  TWICE  ThE  KINETIC  ENERGY 

3  WSB=AMX ( I ) * ( U ( I ) **2+V ( I ) **2+W (I)**2)+ 

1 AMX ( 1+1 >  * ( U ( 1+1 ) **2+V ( I+l ) **2+W ( 1+ 1 ) **2 ) ♦ 
2AMX(J)*(U(J)**2+V(J)**2+W(J)**2)+AMX(J+l)*(U(J+l)**2 
3+V(J+l)**2+W( J+l)**2) 

WSB=WSB+AMX ( K )  *  ( U  ( K ) **2+V ( K ) **2+W (K) **2) +AMX ( K+l ) * 

1 (U(K+1)**2+V(K+1)**2+W(K+1)**2>+AMX(L)*(U(L)**2+V(L)**2 
2+W ( L ) **2 ) + AMX ( L+ 1 ) * ( U  <  L+l ) **2+V ( L+ 1 ) **  2+ W ( L+ 1 ) **2 ) 

C  CALCULATE  THE  NEW  VELOCITIES 

4  U(KK)=(AMX(I)*U(I)+AMX(I+l)*U(I+l)+AMX(J)*U(J)+AMX(J+l)*U(d+l)+ 
1AMX(K)*U(K)+AMX(K+1)*U(K+D+AMXCL)*U(L)+AMX(L+1)*U(L+1)  )/WSA 

V(KK)=(AMX(I)*V(I)+AMX(I+1)*V(H-1>+AMX(J)*V(J)+AMX<J+1>*V(J+1)+ 
1AMX(K)*V(K>+AMX(K+1)*V(K+1)+AMX(L)*V(L)+AMX(L+1)*V(L+1)>/WSA 
W  (KK)  =  ( AMX ( I ) *W ( I ) +AMX ( 1+1 ) *W ( I+l ) +AMX ( J) *W( J)+AMX t J+l ) *W( J+l )  + 
1AMX(K)*W(K)+AMX(K+1)*W(K+1)+AMX(L)*W(L)+AMX(L+1)*w(L+J.))/WSA 
C  CALCULATE  THE  TOTAL  INTERNAL  ENERGY 

WSC=AMX(I)*AIX(I)+AMX(I+1)*AIX(I+1)+AMX(J)*AIX(J)+ 

1AMX(J+1)*AIX(J+1)+AMX(L)*AIX(L)+AMX(L+1)*AIX(L+1)+AMX(K>* 

2AIX(K)+AMX(K+1)*AiX(K+1) 

P(KK)=0. 

C  SET  the  new  MASSES 

AMX ( KK ) “WSA 

CALCULATE  THE  NEW  KINETIC  ENERGY  (ACTUALLY 
TWICE) 

WS=U(KK)**2+V(KK)**2+W(KK)**2 
E=WSC+WSB/2. 

THE  NEW  SPECIFIC  INTERNAL  ENERGY  IS  THE 
TOTAL  LESS  THE  KINETIC 
AIX(KK>=E/WSA-.5*WS 
AMX(J)=Q. 

AMX (J+l ) =0 • 

AMX(K)=0. 

AMX(K+1)=0. 

AMX(L)=0. 

AMX (L+l )=0 • 

AMX(I+l)=0. 

U(J)=Q. 

U(J+1)=0. 

U(K)=Q. 

U(K+1)=Q. 

U(L)=0. 

U(L+1)=0. 

U(I+1)=0. 

V(J)=0. 

V(J+1)=0. 

V(K)rO. 

V(K+1)=0. 

V(L)=0. 

V(L+1)=0. 


VU  +  1)=0. 

W(J)=0. 

W(J+1)=0. 

W(K)=0. 

W (K+l ) =0 • 

W (L) =0 • 

W ( L+l ) =0 • 

W( I+1)=0. 

AIX(I+1)=0. 

A1X(J)=0.  ‘5 

AIX(J+1)=0. 

AIX(K)=?. 

AIX(K+1 ’=0. 

AIX(L)=0. 

AIX(U+1)=0. 

IFdI-1)  380*380*390 

380  IF ( Jj-1 ) 381*381*390 

381  IF(KK-1)7, 7*390 
390  AMX(1)=0. 

U(I)=0. 

V  ( I )  =0 . 
rid)  =0 
P  ( I )  =0 . 

AIX(I)=0. 

GO  TO  7 

C  CAME  HERE  BECAUSE  OF  ZERO  MASS 

6  AMX(KK)=0. 

U(KK)=0. 

V(KK)=0. 

W(KK)=0. 

P(KK)=0. 

AIX(KK)=0. 

P(KK)=0. 

7  KK=KK+IMAX 
21  I=I+2*IMAX 

C  CHANGE  ALL  CELL  DIMENSIONS 
WS=0. 

DO  10  1=1  *  UMAX 
DX ( I ) =2. 0*DX( I ) 

X ( I ) =WS+DX ( I ) 

WS=X(I) 

10  continue 

i i=i imax+i 
WS=X(IIMAX) 

DO  11  I=II*IMAX 
DX(I)=DX(IIMAX) 

X ( I ) =WS+DX( I ) 

WS=XU> 

11  CONTINUE 

ws=o. 

DO  13  J=1*JJMAX 
DY ( J)=2»0*DY ( J) 

Y(J)=WS+DY(J) 

WS=Y(J) 

13  CONTINUE 
JJ=JJMAX+1 
WS=Y(JJMAX) 


o  o 


00  14  J=JJ'JMAX 
DY( J)=DY  ( JJMAX) 

Y(J)=WS+DY(J) 

WS=Y(J) 

14  CONTINUE 
WS=0. 

DO  16  K=1»KKMAX 
DZ(K)=2.0*D2(K) 

ZCOR(K)=WS+DZ(K) 

WS=ZCOR(K) 

16  CONTINUE 

kk=kkmax+i 

WS=ZCOF<(KKMAX) 

DO  17  K=KK»KMAX 
DZ(K)=DZ(KKMAX) 

ZCOR(k)=WS+DZ(K) 

WS=ZCOR(K) 

17  CONTINUE 
KK=KKMAX+1 

DO  30  K=KK»KMAX 
LL=(K-1)*IXMAX 
DO  30  Ul'IMAX 
L=LL+I 

DO  30  J=1»JMAX 
AMX(L)=0. 

U(L)=0. 

V(L)=0. 

W(L)=0.  \ 

AIX(L)=0. 

P(L)=0. 

30  L=L+IMAX 

DO  100  K=1»KKMAX 
LL=(K-1)*IXMAX 

do  loo  i=i »  Umax 
L=LL+IlMAX+l-I 
M=LL+IMAX-I4+1-I 
DO  100  J=1»JJMAX 
AMX(M)=AMX(L) 

U(M)=U(L) 

V (M) =v (L) 

W(M)=W(L) 

AIX(M)=AIX(L) 

P(M)=P(L) 

AMX(L)=0. 

U(L)=0. 

V(L)=0. 

W(L)=0. 

AIX(L)=0. 

P(L)=Q. 

M=M+IMAX 

l=l+imax 

100  CONTINUE 
13=13/? 

NOW  BEGIN  ADDING  ON  MASS  IN  FRONT »  BOTH 
SIDES  AND  ABOVE. 

II=IIMAX-I4 

JJ=I3+1 


c  NOTE*  I4=N0.  OF  ZOl\|ES  TO  THE  RIGHT  TO  ADD. 

C  13=  INITIAL  INTRFACE  BETWEEN  PROJECTILE  AND  TARGET. 

DO  200  K=1»KMAX 
LL=(K-1)*IXMAX 
DO  2u0  1=1*11 
L=LLH  +  I3*IMAX 
DO  200  J=JJ» JMAX 
A MX ( L) =DX ( I > *DY ( J ) *DZ ( K ) *RHONOT 
200  L=L+IMAX 
II=lMAX-I4 
IL= I I MAX- 14+1 
JL= JJMAX+1 
DO  300  K=1»KMAX 
LL=(K-1)*IXMAX 
DO  300  I=IL* I I 
L=LL+I+JJMAX*IMAX 
DO  300  J=JL*JMAX 
AMX(L)=DX(I)*DY(J)*dZ(K)*RHONOT 
'300  L=L+IMAX 

II=IMAX-I4+1 
JJ=I3+1 

DO  400  K=1*KMAX 
LL=(K-1)*IXMAX 
DO  400  1=1 1  * IMAX 
L=LL+I+I3*IMAX 
DO  400  J=JJ» JMAX 
AMX ( L ) =DX ( I ) *DY { J ) *DZ ( K ) *RHONOT 
400  L=L+IMAX 
I I=lMAX-I4 
IL=IIMAX-I4+1 
JJ=I3+1 

IF (JJ-J JMAX) 70 0*70 0*800 
700  KL=KKMAX+1 

DO  500  K=KL*i\MAX 
LL=(K-1)*IXMAX 
DO  500  I=IL* I I 
L=LL+I+I3*IMAX 
DO  500  J=JJ» JMAX 
AMX ( L) =DX ( I ) *DY ( J) *DZ (K ) *RHONOT 
500  L=L+IMAX 
800  WS=T+0TNA 

C  RESET  ACTIVE  GRID  COUNTERS 

11=14-1 
I 2= I MAX- 14+2 
Jl=Jl 

J2=JJMAX+2 
K1=K1 

K2=KKMAX+2 
NK=NC+1 

WRITE (6* 9000 )WS»Nk 
9000  FORMAT ( 1H  // //22H  PROBLEM  REZONED  AT  T=1PE12.6#6X,5HCYCLEI4////) 
WRITE (6*8000) IMAX, (X( I) *1=1 *IMAX) 

WRITE(6»8003) IMAX* (DX ( I ) » 1=1  * IMAX) 

WRITE(6»8001) JMAX, (Y(J) ,J=1»JMAX) 

WRlTt(6, 8004) JMAX, (DY(J) » J=l» JMAX) 

WRITE (6*  8002) KMAX , (ZCOR(K) »K=1*KMAX) 

WRITE (6*8005) KMAX, (DZ(K) #K=1*KMAX) 


000000.000000000000 


73 


WRITE ( 6»  8006) IMAX , JMAX»KMAX , IXMAX »  KMAXA 
RETURN 

Q000  FORMAT ( 1H  /10H  X(I)  1=1 » 12/ ( 5F16 . 6)  ) 

.8001  FORMAT ( 1H  /lOH  Y(j)  J=l» I2/(5F16.6) ) 

8002  FORMAT ( 1H  /13H  ZCOR(K)  K=l> I2/(5F16.6) ) 

8003  FORMAT (1H  /11H  DX(I)  I=l» I2/(5F16.6) ) 

6004  FORMAT ( 1H  /11H  DY(J)  J=l» I2/(5F16.6) ) 

*8005  FORMAT (1H  /11H  DZ (K), K=l» 12/ ( 5F16 . 6) ) 

8006  F0RMAT(7I8) 

END 

01  FOR  ES/S » ES/S » ES/SS 

SUBROUTINE  ES  ES  OOlo 

PH2  0690 


*********  3D0IL  *************************** 

$$$$  FOR  COMPLETE  DETAILS »  SEE  GA-3216» 

METALLIC  EQUATIONS  OF  STATE  FOR  HYPERVELOCITY  IMPACT 
BY  JAMES  TILLOTSON 

if  the  material  is  compressed  »  use  the  condensed 
form  of  the  equation  of  state. 

if  the  material  is  rarefieo  and  if  the  specific 
internal  energy  is  greater  than  e  sub  s  »  USE 
the  rarfied  form  ...  but - 

IF  RAREFIED  and  E  is  LESS  THAN  E  SUB  S» 

USE  THE  CONDENSED  FORM. 

$$SS  NOTE»  NO  NEGATIVE (TENSION)  PRESSURES  ALLOWED 

\ 


10  RHOW=AMX(L)/(DX(I)*DY<J)*DZ(K) ) 
ETA=RHOW/Z ( 33) 

V0W=1 . O/ETA 

11  P1=AIX (L) *RHOW*Z (34) 

12  P2=AIX (L) 

13  P3=Z(3S)*ETA*ETA 

14  P4=Z (36) / (P2/P3+1. 0) *AIX (L) *RHOW 

15  P5=Z(37)*(ETA-1.) 

16  IF(ETA-1.0)50»100»100 
50  IF(V0W-Z(38) )55»55»75 

55  IF(AIX(L)-Z(40) )100»100»75 
75  P7=Z(4l)*(V0w-l.) 

IF(P7-88.0)4002»4002*4003 
4003  P7=88.0 
4002  CONTINUE 
P8=EXP(P7) 

P9=1.0/P8 

P10=2 (42) * ( (VOW-1. )**2) 
IF(P10-88.0)4000»4000»4001 
•  4001  P10=68»0 
4000  continue 

P11=EXP(P10) 

P12=1*0/P11 

P(L)=P1+(P4+P5*P9)*P12 
GO  TO  119 

100  P6=Z(44)*( (ETA-1. )**2) 
P(L)=Pl+P4+P5+P6 


ES  0980 


ES  1010 


ES  1070 


ES  1110 
ES  1120 
ES  1130 
ES  1140 
ES  1150 

ES  117o 
ES  1180 
ES  1190 
ES  1200 
ES  1210 


oooo  ooooo  •  ooooo 


119  IF(P(L) ) 999 » 999 #200 

200  wSGX=.5 

ES 

1260 

GO  To  500 

ES 

1270 

999  P(L)=0. 

wSGX=.5+Z(43) 

GO  To  500 

ES 

1300 

•  bOO  RETURN 

end 

01  FOR  EUIT/S#EUIT/S#EDIT/SS 

subroutine  edit 

EDIT0010 

EDIT0990 

*********  3D0IL  *************************** 

HERE  WE  WILL  DECIDE  WHETHER  TO  HAVE  A  SHORT  PRINT*  LONG  PRINT# 

DUMP  ON  The  binary  tape  or  stop  the  problem 


104  CALL  SLITET(3»K000FX) 

GO  TO(106»108) #K000FX 
106  CALL  SLITE  (3) 

GO  TO  126 

108  IF (CYCLE-CSTOP) 110*122*122 

110  IF ( REZ) 9901 *112*124 

112  IF ( Ay, 0D (CYCLE# DUMPT7 )  )114»124»114 

114  JF(AMOD(CYCLE»PRlNTL) )120»126»120 

120  IF(AMOD(CYCLE»PRlNTS) >140,128»140 

122  CALL  SLITE  (1) 

124  GO  TO  1 

126  CALL  SLITE  (4) 

128  GO  To  6000 
130  GO  TO  1000 
132  CALL  SLITET(4»K000FX) 

GO  T0(134» 136) »KQ00FX 
134  GO  To  5000 


EDIT1040 

EDIT1050 

EDIT1060 

EDIT1070 

EDIT1080 

EDIT1100 

EDIT1150 

EDIT1160 

EDIT1170 

EDIT1180 

EDIT1190 

EDIT1200 

EDIT1210 

EDIT1220 

EDIT1230 


CALCULATE  THE  ENERGY  CHECK  =  (ETH-EJ/ETH  AT  CYCLE  N+l  - 
(ETH-E)/ETH  AT  CYCLE  M  ALL  DIVIDED  THRU  BY  THE  NUMBER  OF  CYCLE 
BETWEEN  ENERGY  CHECKS  =(CYCLE(N+1)  -CYCLE  M) 


136  IF(ABS(ECK)-DMIN) 140#140»9905  EDIT124o 

140  CALL  SLITET(1»K0Q0FX)  EDIT1250 

GO  T0C142*144) »K000FX  EDIT1260 

142  REWIND (N7) 

CALL  SLITE  (1)  EDIT1280 

144  GO  TO  10000  EDIT1290 

1  IF(DUMPT7)30»3»3  EDIT1330 

****  DUMP  ALL  THE  CELL-CENTERED  QUANTITIES  r THE  X,S  AND  Y,S 
AND  ZCOR » S  AND  THE  entire  Z  BLOC’K 

3  bACKSPACE  N7 

wS=5S5.0  EDIT1360 


WRITE (N7)WS»CYCLErN3 
WRITE(N7) (Z(L) »L=1»MZ) 

WRITE (N7) (UCK) »V(K) »W(K) »AMXCK) *AIX(K) »K=1*KMAXA> 

WRITE (N7) ( X ( X ) » I=1»IMAX) 

WRITE (N7) (Y(J) »J=lrUMAX) 

WRITE(N7) (ZCOR(K) »K=1*KMAX> ' 

WS=666.0  EDIT1480 


WRITE (N7)WS»WS#WS 
WRITE (6*  8120) NC 
.  30  60  TO  126 

6000  CONTINUE 

DO  6012  I=l»18 
6012  PR ( I )  =0  • 

C 

C  CALCULATE  THE  TOTAL  INTERNAL  AND  KINETIC  ENERGY  »AND  THE 
C  TOTAL  MASS. 

C 

00  6028  K=1»KMAXA 
6017  PR(1)=0.0 
PR (2) =0.0 
PR (4 ) =0 • 

6019  IF(AMX(K) )9917#6028»6020 

6020  WSB= (U (K ) **2+V (K) **2+W ( K ) **2) * .5 
PR ( 5) =PR ( 5) +AIX (K ) *AMX (K) 

PR(6)=PR(6)+WSB+AMX(K) 

PR(8)=PR(8)+AMX(K) 

6028  CONTINUE 

PR ( 3 ) =PR ( 1 ) +PR ( 2 ) 

PR (7) =PR (5) +PR (6) 

XNRG=PR(7) 

PR ( 9 ) =PR ( 1 ) +PR ( 5 ) 

PR ( 10 ) =PR (2) +PR (6) 

PR ( 1 1 ) =PR ( 3 ) +PR ( 7 ) 

PR(12)=PR(4)+PR(8) 

WSA=(ETH-PR(11)  J/ETH  •' 

IF ( C YCLE ) 9931 » 9931*9932 

9931  NPC=l 

9932  PR(18)=(WSA-DNN)/FL0AT(NPC) 

ECK=PR (18) 

DNN=WSA 

NPC=0 

SUM=O,0 

C  RAOEB=  TOTAL  POSITIVE  Z  MOM. 

C  RADER  =  TOTAL  POSITIVE  X  MOM. 

C  RADE1 =  TOTAL  POSITIVE  Y  MOM. 

SUMBzO. 

SUMR=0. 

SUMTzO  .  '• 

DO  810  K=1»KMAXA 
IF(AMX(K) )  010 » 310 »  802 

802  IF(V(K) )804»804»803 

803  SUMT=SUMT+AMX(K)*V(K) 

804  IF(U(K) )805>805»806 
806  SUMR=SUMR+AMX(K)*U(K) 

805  IF(W(K) )810»810#808 
808  SUMB=SUMB+AMX(K)*W(K) 

810  CONTINUE 

RAUEB=SUMB 

RADERzSUMR 

RADETzSUMT 

WRITE <6#8116)PR0B»NC*TfDTNA#TRAD»DTRAD»NR>Nl»N2»N3»N4 
WRITE (6» 8117) (PR(I) » 1=1*8) 

WRITE (6» 8118) (PR(I) #1=9# 12) 

WR I TE ( 6  »  8 119 ) RAOEB » RADER  »  R ADET  *  UVmAX  * ETH » ECK 


EDIT1510 


ED.T1760 

EDIT1770 

EDIT1790 

EDIT1910 

EDIT1920 

EDIT1930 

EDIT1940 

EDIT1950 

EDIT1960 

EDIT1970 

EDIT1980 

EDIT1990 

EDIT2000 

EDIT2010 

EDIT2020 


EDIT2040 
•  EDIT2050 
EDIT2060 
EDIT2070 


oooooooooooooooooooo  .  ‘  .  ooo 


9042 

WRITEtfcf  9040>Nl0f  NllfN9»U»I2»  Jl*  J2»K1*K2 

WRITE ( b » 9042 ) AMLOSTt ELOST, XMLOST »  YMLOST » ZMLOST t EN£G 

FORMAT l 1P7E14.7) 

6090 

00  To  130 

EDIT2410 

1000 

GO  TO  1030 

EDIT2460 

1030 

WRITe(6>6116)PR0B#NC»T>DTNA»TRAD#DTRAD'NR»N1»N2#N3»N4 

JMAXsjMAX 

EDIT2480 

5000 

WRITE (o» 6307 ) XI » X2» XMAX» Yl » Y2 » Y (JMAX) 
oO  To  132 

WRITE16»B116)PR0B»NC»T»DTNA»TRAD»DTRAD»NR»N1»N2»N3»N4  • 

notice  the  limits  of  the  do  loops 

DO  1126  KK=K1»K2 

HERE  WE  PREPARE  FOR  THE  LONG  PRINT 

WRlTE(6»9041)KKrZC0R(KK) »DZ(KK) 

5004 

DO  5050  1=11 » 12 

CALL  SLITE  (4) 

EDIT2930 

J=J2+1 

K=J2*IMAX+I+ (KK-1 ) +IXMAX 

DO  5046  L=J1 t J2 

J=J-1 

EDIT2970 

K=K-IMAX 

EDIT2980 

5012 

IF ( AMX (K )) 9917 f5046i 5014 

EDIT2990 

5014 

CALL  SLITET(4»K000FX) 

EDIT3000 

GO  T0(5016f 5019) ,K000FX 

*  EDIT3010 

5016 

WRITE(6»Q135) I»X(I) »DX(I) 

5019 

WS=AMX(K)/(DX(I)*DY(J)*DZ(KK)) 

WSC=p(K)*l.E+4 

EDIT3050 

J=ROw  NUMBER 


U=  X  component  of  velocity  cm./sh. 

V=  Y  COMPONENT  OF  VELOCITY  CM./SH. 

WSC  =  PRESSURE  IN  MEGABARS 

AMX  =  MASS  IN  GRAMS 

WS=  DENSITY  IN  GRaMS/CM  CUBED 

W=  Z  COMPONENT  OF  VELOCITY  CM./SH. 

AIX=  SPECIFIC  INTERNAL  ENERGY  IN  JERKS/GM.  (1  JERK=  10<16>ERGS> 
Y=  TOP  COORDINATE  OF  THIS  ROW  IN  CM. 


5016 

WRITE (6» 8108) J#U(K) »V(K) »WSC» AMX (K) » WSr AIX (K) r W (K) * Y ( J) 

5046 

CONTINUE 

EDIT3080 

5050 

continue 

EDIT3090 

1126 

CONTINUE 

* 

GO  To  136 

EDIT3100 

9901 

NK=110 

EDIT3150 

GO  TO  9999 

EDIT3160 

9905 

NK=136 

EDIT3180 

r>  r> 
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GO  TO  9999 
9917  NK=601b 

GO  TO  9999 
'  9999  NR=6 

WRITL(6»8002)I» J»K»KP»I1»I2»NK»NR 
*002  FORMAT (815) 

CALL  UNCLE 
CALL  DUMP 

10000  return 


FORMATS 


EDIT3190 

EDIT321Q 

EDIT3220 

EDIT3320 


EDIT3340 

EDIT3350 

EDIT3360 


0108  FORMAT ( I3» 1X» 1P8E12. 5) 

611to0FORMAT(bHlPR0bL£MLX»5HCYCLE9X#4HTIMEl3X»2HDT13X,4HTRADllX»5HDTRADlEDlT3380 
12X » 2HNR6X » 2HN14X  #  2HU24X  #  2HN34X  #  2HN4/  (F7.1»Ill»2X»  1P4E16 . 7 , 1 10 » 2X  *  4EDIT3390 
216))  EDIT3400 

81 17  OF  OHM  AT  ( 1H0//17X2H/..  I16X»  2RAK14X#  5HAI+AK15X»  2HAM/4H  D0T3X» 1P4E18.7/3EDIT3410 
1H  X4X»4C18.7)  EDIT3420 

81l80FORMAT(12X»l3H - 5X#  13H - - — - — 5X»  l3H - 5EDIT3430 

1X»13h - /7H  T0TALS1P4E18.7)  EDIT3440 

81190FORMAT(2HO  //16X# 5HRADEbl3X#5HRADERl3X» 5HRADET12X»7HMAX  VEL13X»3HTEDIT345o 
1HE12X»9hREL  LRR0R/7X»1P6E18.7////)  EDIT3460 

8120  FORMAT(1HO//21H  TaPE  7  DUMP  ON  CYCLEIb////)  EDIT3470 

81350FORMAT  (lri  ///4H  I  =I3#6X#6HX(  I)  =F12. 3» 6X» 7hDX ( I)  =F12.3//3H  J8X»EDIT3520 
11HX1GX* 1HY10X»  3HF/A9X#  3HAMX9X#  3HRH08X 1 3HAIX9X» 4H  W  8X»2H  Y/) 

8307  F0RMaT(5H  XI  =lPEl2,b# 3X #4HX2  =E12.5» 3X»6HXMAX  =E12.5»6X»4HY1  =E12 
1.5»3X»4HY2  =E12.5#3X»6HYMAX  =E12.5) 

9040  FORMATUH  /  916)  \ 

9041  FORMAT (1M  ///4H  K  =I3»6X»9HZC0R(K)  =F12.3»6X»7MDZ(K)  =F12.3) 

END 


EDIT3630 
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5.  INPUT  AND  DEFINITIONS  OF  THE  VARIABLES 


5*1  Normal  Input  for  the  TRIOIL  Code 

An  asterisk  (*)  implies  that  the  data  on  the  card  is  to  he  converted 
to  fixed  point  data  (requires  a  2  punch  in  Column  l).  All  data  loaded  via 
the  card  routine  is  read  by  a  floating  point  format.  A  double  asterisk  (**) 
signifies  that  this  is  the  last  data  card  in  this  set,  and  requires  a  one 
in  Column  1. 


First  Set 

The  number  of  BCD  (header  cards)  that  will  be  read  in  (Columns  1-3, 
format  13). 

N  BCD  cards,  alphameric  and  or  numeric  in  Columns  2-72. 


Second  Set 

Location 

Symbol 

Description 

*103 

N7 

Binary  tape  number  (data  tape) 

**36271 

PK(1) 

Problem  number 

36272 

PK(2) 

The  cycle  number  to  start  the  calculation 

36273 

PK(3) 

If  <  0.  code  assumes  that  this  is  a  re¬ 
start  or  that  the  CLAM  code  has  generated 
the  initial  data.  (NOTE:  At  this  time  ad 
writing  the  report,  a  three-dimensional 
version  of  the  CLAM  code  is  not  available. 
If  PK(3)  the  code  will  call  subroutine 

set-up. 

Third  Set 

Location 

Symbol 

Description 

1 

PROB 

Problem  number,  identical  to  the  value  in 
PK(1) 

3 

DT 

The  time  step.  Atn  in  shakes.  (NOTE: 

1  shake  =  10"°  seconds) 

4 

PRINTS 

Short  print  frequency  in  cycles 

5 

PRINTL 

Long  print  frequency  in  cycles 

6 

DUMFT7 

Binary  tape  dump  frequency  in  cycles 
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Location 

Symbol 

Description 

✓ 

• 

7 

CSTOP 

Cycle  at  which  the  problem  will  stop 

• 

• 

13 

FFA 

Upper  limit  for  stability  and  to  calculate 

At  if  CABLN  -  0. 

• 

14 

FFB 

Lower  limit  for  stability  and  to  calculate 

At  if  CABLN  -  0. 

20 

EMIN 

If  ECK  (Z(24))  >  DMIN,  problem  will  stop 
because  of  poor  energy  consertation 

-4 

~  10  pQ  =  minimum  density  for  mass  trans¬ 
port  at  a  free  surface  within  the  grid 

23 

TQZONE 

25 

SBOUND 

=  1.0,  fraction  of  A  in*  the  weighted 
velocity  term  in  the  calculation  of  the 
mass  flux 

26 

CABLN 

'if  <  0.  the  code  will  control  At  at 

PCSTAB  (an  input  number)  of  stability 

•  • 

• 

• 

• 

NOTE: 

These  2  options  re- , 
quire  that  you  load> 
At  (location  3)  at 
time  t  ■  0. 

If  =>  0.  the  code  will  control  At,  de¬ 
creasing  At  if  |u  or  c|  At/Ax  or 
|v  or  c[  At/Ay  or  |w  or  c|  At/Az  ex¬ 
ceed  FFA  (an  input  number)  and  increasing 

At  if  the  above  terms  are  less  than  FFB 
(an  input  number) 

If  >  0.,  the  At  that  is  loaded  at  t  -  0., 
will  remain  constant,  regardless  of  any 
^stability  considerations 

29 

WSGD 

=  GAMMA  for  the  dot  material 

30 

WSGX 

■  GAMMA  for  the  x  material 

45 

DTCHK 

.4 

~  10  pQ,  any  cell  that  has  a  density  less 
than  this  value,  will  be  bypassed  for  stab¬ 
ility  checks 

46 

PCSTAB 

~  .25,  fraction  of  stability  as  determined 
by  the  Courant  condition  or  particle  velocity 

58 

2(58) 

Initial  x  velocity  component  of  the  pro¬ 
jectile  in  cm/shake 

59 

2(59) 

Initial  z  velocity  component  of  the  pro¬ 
jectile  in  cm/shake 

66 

RHONOT 

Initial  density  (pQ)  of  all  material  (since 
this  is  a  one-material  code) 

• 

• 

67 

VELOC 

Initial  y  velocity  componenet  of  the  pro¬ 
jectile  in  cm/shake 

• 

68 

BUG 

~  .05,  epsilon  for  determining  whether 

special  features  will  be  used  to  empty  the 
bottom  cells  in  the  projectile 
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Location 

69 

70 

*36 

*87 

*88 

75 

76 

77 

78 

79 

80 

*95 

*96 

*97 

*98 

*99 

*100 

*101 

*102 

*103 

36151 

36181 

36211 

50 

51 

52 

53 


Symbol 


Description 


Z(69) 

Z(70) 

iMAX 

JMAX 

MX 

Z(75) 

Z(76) 

Z(77) 

z(78) 

z(79) 

z(80) 

13 

14 

N1 

N2 

K3 

N4 

N5 

N6 

N7 

Dx 

Dy 

Dz 

51 

52 

53 

54 


=  percent  of  (PCSTAB)  at  time  t  B  0.  vised, 
for  problems  in  which  most  of  the  energy 
at  t  =  0.  is  in  the  form  of  internal  energy 

Factor  to  increase  Z(69)  every  cycle  up 
until  Z(69)  is  equal  to  1.0 

Maximum  number  of  zones  in  the  x  direction 

Maximum  number  of  zones  in  the  y  direction 

Maximum  number  of  zones  in  the  z  direction 

Density  of  material  leaving  the  left  boundary 
of  grid  in  order  to  trigger  rezone 

Similar  term  for  the  front  boundary  of  grid 

Similar  term  for  the  right  boundary  of  grid 

Similar  term  for  the  back  boundary  of  grid 

Similar  tern  for  the  top  boundary  of  grid 

Similar  term  for  the  bottom  boundary  of  grid 

The  original  (j)  value  of  projectile- 
target  interface 

The  number  of  zones  to  add  on  to  the  right 
side  of  grid  after  rezone.  i4  +  the  number 
to  add  on  the  left  is  *  to  iMAX/2 

If  =  0.  the  left  side  of  grid  is  transmittive, 
otherwise  reflective 

Similar  term  for  the  right  side  of  grid 

Similar  term  for  the  top  side  of  grid 

Similar  term  for  the  bottom  side  of  grid 

Similar  term  for  the  back  side  of  grid 

Similar  term  for  the  front  side  of  grid 

Binary  tape  number 

=  Ax  to  be  used  for  all  (i) 

=  Ay  to  be  used  for  all  (j) 

=  Az  to  be  used  for  all  (k) 

=  interface  (j)  value  +  1,  between  the 
projectile  and  the  target 

=  back  (k)  boundary  +  1  of  the  projectile 

=  front  (k)  boundary  of  the  projectile 

=  the  left  (i)  boundary  +  1,  of  the  pro¬ 
jectile 
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Location 

Symbol 

Description 

54 

S5 

=  the  right  (i)  boundary  of  the  projectile 

55 

s6 

=  the  bottom  (j)  boundary  +  1,  of  the  pro¬ 
jectile 

56 

S7 

=  the  top  boundary  of  the  projectile  (j) 

*93 

11  > 

=  minimum  and  maximum  values  of  the  active 

*94 

12  J 

grid  in  the  x  direction 

*108 

k! 

=  minimum  and  maximum  values  of  the  active 

*109 

k2  J 

grid  in  the  z  direction 

*110 

J1  ^ 

=  minimum  and  maximum  values  of  the  active 

*111 

J2  J 

grid  in  the  y  direction 

47 

Z(47) 

■  A°  c  *  c  +  A  (P)*2 

48 

Z(48) 

49 

Z(49) 

-A  j  °  1 

Where  P  is  in  megabars  and  cQ  and  A-^ 
are  in  units  of  10'  cm/sec,  converted  to 
cm/shake  in  CDT  routine 


33 

z(33) 

«  Po  1 

34 

Z(34) 

-  a 

35 

z(  35 ) 

=  Eo 

36 

Z(36) 

■  b 

37 

Z(37) 

=  A 

38 

z(38) 

-  Vs 

v  Constant  for  metallic  equation 

39 

z(39) 

Blank 

*  of  state  (Tillotson  formulation) 

40 

z(4o) 

-  Es 

4l 

z(4i) 

=  a 

42 

Z(42) 

-  0 

43 

Z(43) 

Blank 

44 

2(44) 

=  B  J 

16 

xMAX 

=  epsilonics  on  the  velocity,  if  |  u|  or 
| v|  or  |w|  <  e,  it  is  set  to  0.  and  the 
books  are  kept 

15 

TMASS 

=  epsilonics  on  the  specific  internal  energy 
If  I  <  e,  it  is  set  to  0.  and  the  books 
are  kept 

**22 

REZFCT 

No  meaning  in  this  code 

Last  Set 

Location 

Symbol 

Description 

**22 

REZFCT 

No  meaning 

« 


1 


t 


I 
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5.2  List 

of  Common 

for  TRIOIL 

Symbol 

Location 

No.  of 
Words 

Units 

Description 

AIX 

151 

6000 

jerks/g 

Specific  internal  energy  for 
cell  (L) 

AMX 

6151 

6000 

g 

Mass  in  cell  (L) 

U 

12151 

6000 

cm/ shake 

X  component  of  velocity  in 
cell  (L) 

V 

18151 

6000 

cm/shake 

Y  component  of  velocity  in 
cell  (L) 

w 

24151 

6000 

cm/ shake 

Z  component  of  velocity  in 
cell  (L) 

p 

30151 

6000 

jerks/cnr 

Material  pressure  in  cell  (L) 

DX 

36151 

30 

cm 

DX(i )  *  X(i )  -  X(i-l) 

BY 

36181 

30 

cm 

DY(J)  =  Y(j)  -  Y(j-l) 

BZ 

36211 

30 

cm 

DZ(k)  =  ZCOR(k)  -  ZCOR(k-l) 

UL 

36241 

30 

cm/shake 

Velocity  at  the  left  of  cell 
in  PHI 

FLEFT 

36241 

30 

g  cm/shake 

X  momenta  of  mass  crossing 
left  side  of  cell  (PH2) 

PL 

36271 

30 

O 

jerks/cnr 

Temporary  pressure  array  at 
the  left  interface  in  PHI 

PK 

36271 

30 

none 

Not  used  (except  input) 

YAMC 

36271 

30 

g  cm/shake 

Y  momenta  of  mass  crossing 
left  side  of  cell  PH2 

X 

36301 

30 

cm 

X(i)  *  right  dimension  of 
cell  (L) 

Y 

36331 

30 

cm 

Y(j)  =  top  dimension  of  cell 
(L) 

ZCOR 

36361 

30 

cm 

ZCOR(k)  »  Z  or  front  dimension 
of  cell  (L) 

PR 

36391 

50 

many 

Used  for  editing  in  the  EDIT 
routine 

SIGC 

36441 

30 

jerks/g 

=  specific  energy  of  the  mass 
crossing  left  side  of  cell  (PH2) 

GAMC 

36471 

30 

g 

Mass  crossing  left  side  of 

cell  in  (PH2) 


No.  of 


Symbol 

Location 

Words 

Units 

Description 

ZMOM 

36501 

30 

g  cm/shake 

Z  momenta  of  mass  crossing 
left  isde  of  cell  (PH2) 

BXMOM 

36531 

700 

g  cm/shake 

X  momenta  of  mass  crossing 
back  side  of  cell  (PH2) 

UBIND 

36531 

700 

cm/ shake 

Velocity  at  back  interface 
of  cell  in  PHI 

PBIND 

37231 

TOO 

Jerks/cm^ 

Pressure  array  at  back  in¬ 
terface  of  cell  in  (PHl) 

BMASS 

37231 

700 

g 

Mass  crossing  back  interface 
of  cell  in  PH2 

BYMOM 

37931 

700 

g  cm/shake 

Y  momenta  of  mass  crossing 
back  surface  in  PH2 

BZMQM 

38631 

700 

g  cm/shake 

Z  momenta  of  mass  crossing 
back  surface  in  PH2 

BENR 

39331 

700 

Jerks /g 

Specific  energy  of  this  mass 
crossing  the  back  surface 
in  PH2 

AREA 

40031 

1 

none 

Flag  in  PH2  for  the  Z 
direction 

BIG 

40032 

1 

none 

Not  used 

BOUNCE 

40033 

1 

none 

Not  used 

PABOVE 

40034 

1 

Jerks/cnr 

Pressure  at  the  top  of  cell 
(L)  in  PHI 

PBLO 

40035 

1 

Jerks/cnr 

Pressure  at  the  bottom  of 
cell  (L)  in  PHI 

PiDTS 

40036 

1 

none 

Not  used 

PRR 

40037 

1 

Jerks/cnr 

Pressure  at  the  right  of  cell 
(L)  in  PHI 

RHO 

40038 

1 

none 

Not  used 

SIG 

40039 

1 

cm 

Minimum  Ax,  Ay  or  Az  for  a 
cell  in  CDT  routine 

UVMAX 

40040 

1 

shake  ^ 

(Maximum  velocity )/minimum 
(Ax  or  Ay  or  Az) 

VABOVE 

40041 

1 

cm/shake 

Velocity  at  the  top  of  cell 
(L)  in  PHI 

VBLO 

40042 

1 

cm/shake 

Velocity  at  the  bottom  of 
cell  (L)  in  PHI 

VEL 

40043 

1 

none 

Maximum  (a-l)  in  the  CDT 
routine,  a  flag  for  subcycling 
in  PHI  and  a  flag  for  the  y 
direction  in  PH2 
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No.  of 


Symbol 

Location 

Words 

Units 

Description 

WPS 

40044 

1 

none 

Working  storage 

WS 

40045 

1 

none 

Working  storage 

USA 

40046 

1 

none 

Working  storage 

WS3 

4004? 

1 

none 

Working  storage 

use 

40048 

1 

none 

Working  storage 

3_ 

40049 

1 

none 

Temporary  indices 

ii 

40050 

1 

none 

Temporary  indices 

IN 

40051 

1 

none 

Temporary  indices 

iR 

40052 

1 

none 

Temporary  indices 

iWS 

40053 

J. 

none 

Temporary  indices 

iWSA 

40054 

1 

none 

Temporary  indices 

iWSB 

40055 

1 

none 

Temporary  indices 

iWSC 

40056 

1  . 

none 

Temporary  indices 

J 

40057 

1 

none 

Temporary  indices 

JK 

40058 

1 

none 

Temporary  indices 

jP 

40059 

1 

none 

Temporary  indices 

JR 

40060 

1 

none 

Temporary  indices 

K 

40061 

1 

none 

Temporary  indices 

KDT 

40062 

1 

none 

Flag  in  CDT  for  At  change 

KN 

40063 

1 

none 

Temporary  indices 

KP 

4C064 

1 

none 

Temporary  indices 

KR 

40065 

1 

none 

Temporary  indices 

LRM 

40066 

1 

none 

Temporary  indices 

L 

40067 

1 

none 

Temporary  indices 

K 

40068 

1 

none 

Temporary  .....dices 

MA 

40069 

1 

none 

Temporary  indices 

M3 

40070 

1 

none 

Temporary  indices 

MC 

40071 

1 

none 

Temporary  indices 

MD 

40072 

1 

none 

Temporary  indices 

ME 

40073 

1 

none 

Temporary  indices 

MZ 

40074 

1 

none 

Set  -  150  in  input  required 
in  EDIT 

N 

40075 

1 

none 

Temporary  indices 

Symbol 

No.  of 

Location  Words 

Units 

Description 

REZ 

40076 

1 

none 

Not  used 

TRAD 

40077 

1 

none 

Not  used 

DTRAD 

40078 

1 

none 

Not  used 

RADEB 

40079 

1 

none 

Total  positive  Z  momentum 

RADER 

40030 

1 

none 

Total  positive  X  momentum 

RADET 

40081 

1 

none 

Total  positive  Y  momentum 

XI 

40082 

1 

none 

Not  used 

X2 

40083 

1 

none 

Not  used 

Y1 

40084 

1 

none 

Not  used 

Y2 

40085 

.1 

none 

Not  used 

iMAXA 

40086 

1 

none 

Not  used 

The  Z  Block 

Location 

Symbol 

Units 

Description 

Z(l) 

PROB 

none 

Problem  number 

Z(2) 

CYCLE 

none 

Floating  point  value  of  the 
cycle  number 

z(3) 

DT 

shake 

&t  hydro  =  tn  -  tn  ^ 

Z(4) 

PRINTS 

none 

Cycle  frequency  for  short  print 

Z(5) 

PRINTL 

none 

Cycle  frequency  for  long  print 

z(6) 

DUMPT7 

none 

Cycle  frequency  for  binary 
tape  dump 

Z(7) 

CSTOP 

none 

Cycle  number  at  which  problem 
stops 

Z(8) 

PIDY 

none 

-  it  -  3-1415927 

Z(9) 

GAM 

none 

Not  used 

z(io) 

GAMD 

none 

=  l./(or-l.)  computer  in  input 

Z(ll) 

GAMX 

none 

=  l./(ax-l. )  routine 

Z(12) 

ETH 

Jerks 

=  total  energy  in  the  system 
(originally  set  =  to  (iMAX) 

( JMAX ) ( kMAX ) 

£1  ™(l>[x<l)Hj^l)"(l>] 

Changed  to  PHI  at  transmittive 
boundaries  and  in  PH2  if  mass 
leaves  the  system 
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Location 

Symbol 

Units 

Description 

z(l3) 

FFA 

none 

Upper  limit  for  stability  and 
used  to  calculate  At  only  if 
CABLN  =  0. 

Z(14) 

FFB 

none 

Lower  limit  for  stability  and 
used  to  calculate  At  only  if 
CABLN  =  0. 

Z(15) 

mss 

none 

Epsilonics  on  minimum  specific 
internal  energy 

z(i6) 

XMAX 

none 

Epsilonics  on  minimum  velocity 
components 

z(17) 

MX 

none 

Not  used 

z(l8) 

ZMAX 

none 

Not  used 

z(  19) 

DM 

none 

=  (ETH-E^'^/ETH)  for  energy 
check 

Z(20) 

DMIN 

none 

if  ECK  (see  definition  in  Z(24) 
is  >  DMIN,  problem  will  stop 
because  of  energy  violation 

Z(21) 

DTNA 

shake 

Atn_1 

Z(22) 

REZFCT 

none 

Not  used 

z(23} 

TOZONE 

g/cm3 

If  the  mass  flow  across  a  free 
surface  within  the  grid  produces 
a  density  <  TOZONE,  the  mass 
flow  is  set  to  zero 

Z(24) 

ECK 

none 

Is  the  energy  check  (ETH-En/ETH) 

-  (eth-en_kpc/eth)  /npc 

Where  NPC  =  cycle  frequency  at 
which  the  energy  check  is  made 

Z(25) 

SBOUND 

none 

Fractions  of  A  in  mass  weighted 
velocity  (suggested  number  *  1.0) 

Z(26) 

CABLN 

none 

Defined  in  Section  5*1 

Z(27) 

T 

shake 

Total  time  up  to  cycle  NC 
tn+1  =  tn  +  At 

z(28) 

GMAX 

none 

Maximum  of  the  two  gammas 
(7X  or  7.) 

z(29) 

WSGD 

none 

7. 

z(30) 

WSGX 

none 

7X  and  (7j*iax_i,  )  in  CDT  routine 

z(3D 

GMADR 

none 

%/(7-l* ) 

z(32) 

GMAXP 

none 

7X/(7X-1-) 

m 


Location 

Z(33) 
Z(3*0 
Z  ( 35 ) 
Z(36) 
Z(3T  j 
Z(38) 
Z(39) 

Z(UO 

Z(4l) 

Z(42) 

Z(43) 

Z(U4) 

Z(45) 


Z(46) 

Z(4T) 

Z(U8) 

Z(49) 

z(50) 
Z(51) 
Z(52) 
Z(  53 ) 
Z(5M 
Z(55) 
Z(56) 

Z(57) 

Z(58) 

Z(59) 
Z(60) 
Z(6x) 
Z(62) 
Z(63) 
Z(64) 
Z(65) 
Z  (66) 

Z(6t) 

z(68) 
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Symbol 

Units 

Description 

p„ 

g/cm^ 

o 

a 

none 

Eo 

jerks/g 

b 

A 

V. 

none  , 
Jerks/cnr 
none 

For  metallic  equations  of  state 
>  (Tillotson  formulation) 

s 

none 

none 

Es 

Jerks/g 

a 

none 

0 

none 

0 

none  , 

B 

Jerks/cnr  J 

DTCHK 

g/cnr 

Density  checks,  if  P(l)  < 
the  stability  check  for  cell  (L; 

will  be  bypassed 

PC STAB 

none 

$  of  stability,  used  if  CABLN  is 
<0.  A  recommended  value  is  .25 

CNOT 

BFACT 

10"3  cm/sec 
none 

\  c  *  speed  of  sound  =  CNOT  + 

/  BFACT (P(l))EPS  where  P(L)  is  in 

EPSi 

none 

m 

J  megabars 

SI 

none  " 

\ 

S2 

none 

S3 

none 

V  Used  in  set-up  routine 

S4 

none 

/ 

S5 

none 

s6 

none 

S7 

none 

s8 

none 

Not  used 

S9 

none 

Initial  X  velocity  component 
of  the  projectile  in  cm/shake 

sio 

none 

Similar  term  for  the  Z  direction 

AMLOST 

g 

Mass  thrown  away  (PH2) 

ELOST 

jerks 

Energy  thrown  away  (PH2) 

XMLOST 

g  cm/shake  Total  X  momenta  thrown  away  (PH2) 

YMLOST 

g  cm/shake  Total  Y  momenta  thrown  away  (PH2) 

ZMLOST 

g  cm/shake  Total  Z  momenta  thrown  away  (PH2) 

ENEG 

jerks 

Energy  added  to  system  if  I  <  0. . 

RHONOT 

VELOC 

g/ cnf 
cm/ shake 

Initial  density  of  material 

Initial  velocity  of  pellet  in 
the  Y  direction  (cm/shake) 

BUG 

none 

Epsilon  for  emptying  pellet  (~  .01) 
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Location 

Symbol 

Units 

z(69) 

none 

Z(70) 

none 

Z(71) 

none 

Z(72) 

none 

Z(73) 

none 

Z(7t) 

none 

Z(75) 

none 

z(76) 

none 

Z(  77 ) 

none 

Z(78) 

none 

z(79) 

none 

Z(80) 

none 

Z(8l) 

NPR 

none 

Z(82) 

NPRi 

none 

Z  (83 ) 

NC 

none 

Z(84) 

NPC 

none 

Z(85) 

NRC 

none 

z(86) 

iMAX 

none 

Z(87) 

jMAX 

none 

Z(88) 

kMAX 

none 

Z(89) 

kMAXA 

none 

z(90) 

iXMAX 

none 

z(9i) 

NOD 

none 

z(92) 

NOPR 

none 

z(93) 

il 

none 

z(9M 

i2 

none 

z(95) 

13 

none 

Description 

Defined  in  Section  5*1 
Defined  in  Section  5*1 
j  value  of  top  of  target 
Not  used 
Not  used 
Not  used 

Defined  in  Section  5*1 

Defined  in  Section  5*1 

Defined  in  Section  5*1 

Defined  in  Section  5*1 

Defined  in  Section  5*1 

Defined  in  Section  5*1 

Index  (working  storage) 

Index  (working  storage) 

Cycle  number  in  fixed  point 

Number  or  cycles  between  energy 
checks 

Index 

Maximum  number  of  zones  in  X 
direction 

Maximum  number  of  zones  in  Y 
direction 

Maximum  number  of  zones  in  Z 
direction 

Total  number  of  zones  =  (iMAX) 

( jMAX)(kMAX) 

=  ( iMAX)( JMAX) 

Index  (workj  ng  storage ) 

Index  (working  storage) 

Minimum  value  of  i  in  do  loop 
on  i 

Maximum  value  of  i  in  do  loop 
on  i 

Original  interface  between  pro¬ 
jectile  and  target 
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Location 

Symbol 

Units 

Descript?  on 

z(96) 

i4 

none 

Number  of  zones  to  the  right  for 
re zone 

z(97) 

N1 

none 

Flag  at  left,  if  0.,  boundary  is 
transmittive,  otherwise  reflective 

Z(98) 

N2 

none 

Flag  at  right,  if  0.,  boundary  is 
transmittive,  otherwise  reflective 

Z(99) 

N3 

none 

Flag  at  top,  if  0.,  boundary  is 
transmittive,  otherwise  reflective 

Z(100) 

N4 

none 

Flag  at  bottom,  if  0.,  boundary  is 
transmittive,  otherwise  reflective 

Z(101) 

N5 

none 

Flag  behind,  if  0.,  boundary  is 
transmittive,  otherwise  reflective 

Z ( 102 ) 

N6 

none 

Flag  in  front,  if  0.,  boundary  is 
transmittive,  otherwise  reflective 

Z(103) 

N7 

none 

Binary  tape  number  designation 

Z(104) 

N8 

none 

Not  used 

Z(105) 

N9 

none 

k  value  of  zone  that  is  con¬ 
trolling  At 

Z(106) 

N10 

none 

i  value  of  zone  that  is  con¬ 
trolling  At 

z(107) 

Nil 

none 

J  value  of  zone  that  is  con¬ 
trolling  At 

z(lo8) 

kl 

none 

Minimum  value  of  k  in  do  loop 
on  k 

Z(109) 

k2 

none 

Maximum  value  of  k  in  do  loop 
on  k 

Z(110) 

J1 

none 

Minimum  value  of  J  in  do  loop 
on  J 

Z(lll) 

J2 

none 

Maximum  value  of  J  in  do  loop 

on  J 


Z(112)  through  Z(150) 


Not  used 
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